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Abstract 


This  report  describes  NavyFOAM  VI .0,  a  computational  fluid  dynamics  (CFD) 
capability  based  on  Reynolds-averaged  Navier-Stokes  equations  (RANSE)  aimed  at 
predicting  turbulent  single-  and  two-phase  flows  around  ship  hulls.  The  CFD  capability 
employs  a  Finite-volume  discretization  that  allows  use  of  arbitrary  polyhedral  elements.  The 
free  surface  is  captured  using  a  volume-fraction  method  capable  of  accurately  resolving  sharp 
interfaces.  NavyFOAM  has  been  developed  using  an  open-source  CFD  software  tool-kit 
(OpenFOAM)  that  draws  heavily  upon  object-oriented  programming.  The  numerical  methods 
and  the  physical  models  in  the  original  version  of  OpenFOAM  have  been  upgraded  in  an 
effort  to  improve  accuracy  and  robustness  of  numerical  solutions.  The  details  of  NavyFOAM 
V1.0  including  the  numerical  methods  and  the  physical  models  are  described  in  this  report. 
NavyFOAM  V1.0  is  demonstrated  for  a  number  of  flows  including:  underwater  bodies, 
turbulent  free  surface  flows  around  the  DTMB  5415  model  and  the  KVLCC2  double-model. 

It  is  shown  that  the  RANSE  based  approach  can  predict,  with  good  accuracy,  most  of  the 
salient  features  of  the  turbulent  free-surface  flows  around  the  subject  hulls  including 
resistance,  wave  elevation,  hull  boundary  layer  and  wake. 

Administrative  Information 

The  work  described  in  this  report  was  performed  by  the  Computational 
Hydromechanics  Division  (Code  5700)  of  the  Hydromechanics  Department  at  the  Naval 
Surface  Warfare  Center,  Carderock  Division  (NSWCCD).  This  effort  has  been  funded  by  the 
Department  of  Defense  High  Performance  Computing  Modernization  Program  (HPCMP) 
under  the  Computational  Research  and  Engineering  Acquisition  Tools  and  Environments 
(CREATE)  Ship’s  Hydrodynamics  Project. 

Introduction 

In  the  past  two  decades,  computational  fluid  dynamics  (CFD)  has  been  established  as 
an  indispensible  tool  for  design  and  analysis  in  ship  hydrodynamics.  CFD  has  also 
significantly  expanded  its  realm,  covering  a  broad  spectrum  of  applications  including: 
resistance,  powenng,  propulsion,  maneuvering  and  seakeeping.  The  geometrical,  physical 
and  operational  complexity  involved  in  ship  hydrodynamics  applications  has  led  to  the 
addition  of  many  features  and  functionalities  in  CFD  codes.  Furthermore,  CFD  is  frequently 
called  upon  to  tackle  multi-disciplinary  applications  such  as  fluid-structure  interaction  and 
hydroacoustics  applications  that  require  coupling  of  CFD  codes  with  other  computational 
mechanics  software.  Thus,  general-purpose  CFD  codes,  in  attempts  to  cater  to  these  diverse 
needs,  have  become  increasingly  larger  and  more  complex.  Software  complexity  is  a  serious 
issue  which  many  legacy  CFD  codes  face  today,  negatively  impacting  their  overall  efficacy  in 
terms  of  quality  assurance,  packaging,  maintenance  and  extensions. 

The  Department  of  Defense  High  Performance  Computing  Modernization  Program 
(HPCMP)  office,  under  the  CREATE  Ship’s  Hydrodynamics  Project,  has  initiated  an  effort  to 
develop  a  CFD  capability  aimed  at  high-fidelity',  high-performance,  predictions  of 
hydrodynamic  phenomena  occumng  around  surface  ships  and  submarines.  The  ultimate  goal 
of  the  project  is  to  develop  a  high-fidelity  CFD  capability  that  can  drastically  shorten  the  design 


cycles  of  surface  ships  and  submarines,  by  answering  teehmeal  questions  on  various  aspects  of 
hydrodynamic  performance  of  naval  vessels  at  early  design  stages. 

To  meet  the  top-level  requirements  of  the  program,  it  was  considered  imperative  that  the 
new  CFD  software  be  developed  using  modem  software  engineering  practices.  Among  others, 
it  was  concluded  that  object-oriented  programming  (OOP)  with  properly  designed  data 
structure  and  code  architecture  is  essential  to  facilitate  development,  quality  assurance  (QA), 
deployment  (paekaging/release),  maintenance,  and  extension  of  the  software.  Thus,  we  started 
with  OpenFOAM  (Weller  et  al.5),  an  open-source  CFD  software  tool-kit  written  in  C++ 
drawing  heavily  upon  object-oriented  programming  (OOP).  Efforts  to  develop  a  computational 
framework  using  OpenFOAM  had  started  out  earlier  with  propulsors  the  target  applications 
(Kim  et  ah").  The  CREATE  efforts  have  greatly  benefited  from  our  earlier  works  on  turbulence 
modeling,  discretization  schemes  and  solution  algorithms.  As  of  today,  the  OpenFOAM-based 
computational  framework  comprises  a  suite  of  modified  and  newly  written  application  (top- 
level)  solvers  for  single-  and  multi-phase  flows,  utilities  and  physics  libraries  built  around  the 
OpenFOAM  CFD  tool-kit.  We  loosely  refer  to  the  computational  framework  as  “NavyFOAM” 
in  order  to  distinguish  it  from  the  standard  OpenFOAM  offering. 

NavyFOAM  includes  several  top-level  solvers.  Table  1 ,  aimed  at  ship  hydrodynamics 
applications,  sRansFOAM  (single-phase,  steady  RANSE  solver),  ransFSFOAM  (RANSE- 
based  free-surface  solver),  and  ransFSDyMFOAM  (RANSE-based  free-surfaee  solver  with 
moving/deforming  mesh),  to  name  a  few.  That  one  has  to  deal  with  a  number  of  top-level 
solvers  for  different  applications  often  surprises  those  who  are  used  to  the  idea  of  developing 
a  monolithic  CFD  solver  that  can  do  everything.  The  philosophy  adopted  in  OpenFOAM 
eschews  the  monolithic  approach. 

This  report  consists  of  a  number  of  sections,  including: 

•  Technical  description 

•  User’s  Guide 

•  Example  Results 

•  Utility  Programs 

•  Tutorials 

Technical  Description  gives  an  overview  of  the  theoretical  formulation  and  the 
numerical  methods  used  in  the  RANSE  solvers  in  NavyFOAM. 

User’s  Guide  is  intended  to  help  users  learn  how  to  run  the  codes  without  delving  into 
the  details  of  the  implementations.  This  chapter  should  be  considered  as  an  annex  to 
OpenFOAM’s  User’s  Guide.  Those  who  are  interested  only  in  running  the  top-level  solvers 
provided  in  NavyFOAM  should  read  this  chapter  and  the  Tutorials  and  ean  skip  the  other 
chapters  if  they  want  to. 

Example  Results  presents  example  problems  run  with  NavyFOAM  selected  from 
various  applications  including  surface  ships  and  underwater  bodies. 

Utility  Programs  is  an  appendix  that  describes  the  top-level  applications  newly  added 
to  facilitate  post-processing  of  the  CFD  results  obtained  using  NavyFOAM. 
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Tutorials  given  in  the  appendices  provide  step-by-step  instructions  starting  from 
setting  up  the  ease  to  running  the  NavyFOAM  solvers  to  exporting  the  results  for  post¬ 
processing. 


Table  1 .  Examples  of  top-level  RANS  solvers  built  using  the  OpenFOAM  toolkit  for  marine 
propulsor  applications  (GGI:  grid-to-grid  interpolation) 


Solver 

Features/Functionalities 

Applications 

sRansFoam 

Single-phase,  steady,  RANSE,  flow  solver  in  the 
inertial  frame 

Underwater  bodies  (without 
free-surface  effects) 

ransFSFoam 

Two-phase,  unsteady,  RANSE,  flow  solver  in  the 
inertial  or  rotating  frame  with  GGI 

Surface  ships  with  fixed 
sinkage  and  trim 

Propellers  in  open  water 
with  uniform  inflow 

ransFSDyMFoar 

Two-phase,  unsteady  RANSE,  flow  solver  in  the 
inertial  frame  with  dynamic  mesh  motion  with  GGI 

Surface  ships  with 
dynamic  sinkage  and  trim 
prediction 

Propellers  with  non- 
uniform  inflow 
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Technical  Description 

NavyFOAM  employs  a  cell-centered  finite-volume  method  based  on  a  multi-dimensional 
linear  reconstruction  scheme  that  permits  use  of  arbitrary  polyhedral  elements  including 
quadrilateral,  hexahedral,  triangular,  tetrahedral,  pyramidal,  prismatic,  and  hybrid  meshes.  The 
solution  gradients  at  cell  centers  can  be  evaluated  by  applying  the  Green-Gauss  theorem  or  by 
the  least-square  method.  Spatial  and  temporal  discretizations  formally  have  up  to  second-order 
accuracy.  The  volume-fraction  equation  is  solved  using  an  implicit  solver.  The  discretized 
governing  equations  can  be  solved  using  a  choice  of  iterative  linear  solvers  such  as  point-implicit 
Gauss-Seidel  or  algebraic  multi-grid  (AMG)  methods.  Velocity  coupling  to  ensure  mass 
conservation  (continuity)  is  effected  using  a  projection  algorithm.  The  entire  NavyFOAM  solver 
suite  can  be  run  in  parallel  using  domain  decomposition  and  a  public  version  of  MP1  (OpenMPI) 
for  message  passing. 

Governing  Equations 

The  governing  equations  adopted  in  NavyFOAM  consist  of  the  continuity  (mass 
conservation)  equation,  momentum  equations,  turbulent  transport  equations,  and  a  volume- 
fraction  equation.  Which  equations  are  solved  in  a  top-level  solver  depends  on  whether  the  flow 
is  single-phase  or  multiphase. 

Single  Phase  Flow 

For  single  phase  incompressible  flow,  the  governing  equations  consist  of  the  continuity 
equation,  the  momentum  equation,  and  the  turbulence  transport  cquation(s).  The  continuity 
equation  can  be  written  in  a  differential  form  as: 

V-  F  =  0  (1) 

The  momentum  equation  can  be  written  as: 

^-  +  V-(VV)  =  -VP  +  V-{ve/f(vv  +  VpT)} 

where  V  is  the  velocity  vector,  P  =  —  is  the  modified 

P 

p  is  the  density,  veff  a  v  +  v9  is  the  effective  viscosity, 
turbulent  eddy  viscosity. 

Multiphase  Flow 

In  the  volume  of  fluid  (VOF)  method,  the  governing  equations  for  two-phase  flow  consist 
of  the  continuity  equation,  the  momentum  equation,  the  convection  equation  for  volume  fraction, 
and  the  turbulence  transport  equation(s).  The  continuity  equation  is  given  by 

V-P  =  0  (3) 

The  momentum  equation  is  given  by 


(2) 

pressure,  p  is  the  hydrodynamic  pressure, 
v  is  the  kinematic  viscosity,  and  vt  is  the 
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( +  v  •  (pVV)  =  -Vp  +  V  •  {^(vp  +  VpT )}+  pg  +  CTK-V/ 


(4) 


where  ?  is  the  volume  fraction,  g  is  the  gravitational  acceleration  vector,  <X  is  the  surface 
tension  coefficient,  and  K  is  the  interface  curvature,  peJf  =  p  +  //,  is  the  effective  viscosity,  p  is 

the  dynamic  viscosity,  and  pt  is  the  turbulent  eddy  viscosity.  The  density  is  calculated  by 
P  =  YP\  +  (1 -  y)P: .  and  the  dynamic  viscosity  by  p  =  //(//, ,  p: ,  y) .  The  subsenpts  “ l  ”  and  “2” 
refer  to  the  two  phases  or  fluids.  The  convective  transport  equation  for  volume-fraction  is 

§7  +  V-(»V)  =  0  (5) 

Turbulence  Models 

NavyFOAM  allows  users  to  choose  from  the  entire  turbulence  model  suite  available  in 
OpenFOAM.  NavyFOAM  additionally  offers  a  Wilcox’s  k-co  turbulence  model  (Wilcox^),  a 
modified  SST  k-co  model,  and  a  custom  version  of  Spalart  and  Allmaras’  one-equation  model. 
Wall  models  are  implemented  in  these  newly  available  turbulence  models  so  that  the  models  can 
be  used  with  either  a  wall-resolving  (y+  <  1)  or  a  wall-skipping  (y+  >30)  mesh. 

Spatial  and  Temporal  Discretization 

Gradient  Schemes 

Gauss  Integration .  The  gradient  can  be  calculated  using  the  Gauss  theorem 

jv<p</Q=  ^TupdV 

n,  r. 

Assuming  the  gradient  is  constant  in  a  cell,  (6)  can  be  approximated  as 

(7) 

where  the  subscript  P  denotes  the  cell  center,  |  Qe  |  is  the  volume  of  the  cell,  and  Sf  is  the  area 
vector  of  each  face  of  the  cell. 

Cell-Based  Calculation.  In  a  cell-based  approach,  the  face-centcr  value  ^ f  in  (7)  is 
calculated  using  the  cell-center  value  of  ^  in  neighboring  cells,  as  shown  in  Figure  1.  This 
approach  is  used  in  OpenFOAM. 
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Figure  1.  Stencil  in  cell-based  gradient  calculation  in  two-dimensions 

Node-Based  Calculation  In  a  node-based  approach,  the  nodal  value  of  tp  (red  circles  in 
Figure  2)  is  First  calculated  using  all  neighboring  cells  of  the  node,  then  the  facc-ccntcr  value  <pf 

in  (7)  is  calculated  using  nodal  values.  As  a  result,  the  stencil  involved  in  the  gradient  calculation 
(all  blue  squares  in  Figure  2)  is  much  larger  than  that  in  the  cell-based  approach. 


Figure  2.  Stencil  in  node-based  gradient  calculation  in  two-dimensions 

Two  types  of  node-based  gradient  calculations,  i.e.  the  Pseudo-Laplacian- Weighted- 
Averaging  (PLWA)  and  the  Inverse-Distance- Weighted- A veraging  (IDWA)  have  been 
implemented  in  NavyFOAM.  They  differ  in  the  way  that  the  cell-centered  volume  field  is 
interpolated  to  the  node  point  field  More  details  are  described  later. 

Least  Squares.  The  concept  of  least-squares  calculation  of  gradients  is  easily  illustrated 
in  two-dimensions.  There  should  be  no  difficulty  to  extend  it  to  three-dimensions.  Suppose  wc 
want  to  calculate  the  gradient  of  (p  at  the  center  of  cell  /,  see  Figure  3,  the  neighboring  cells  are  k 
=1,  2,  and  3. 
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Figure  3.  Stencil  in  least-squares  gradient  calculation  in  two-dimensions 

In  a  general  form,  let  Nb  be  the  total  number  of  neighboring  cells,  the  neighboring  cell- 
center  value  of  cp  can  be  written  as  a  Taylor  expansion  about  the  center  of  cell  /. 


<Pk  =  9,  +  (v^),  '  A',*  +  Eik  for  /t  =  1 ,  . 


Nh 


(8) 


where  sik  represents  the  higher-order  errors.  Defining  a  total  error  as  the  sum  of  weighted  errors 
using 


E  =  14 


(9) 


A  =1 


where  wik  is  the  weighting  factor.  Omitting  the  index  i  for  brevity,  Equation  (8)  can  be  written  as 

fort=1 . N> 

or  in  the  form  of  Cartesian  components 


Ax* 

1 

s* 

1 

s. 

II 

*■< 

CO 

<Py 

* 

_<P: 

_Az*_ 

(10) 


for  k=  1,...,  Nh 


Substituting  (10)  into  (9)  and  setting 

^  =  0,  «..a  ml  ^-0 

d<px  d<py  d(p. 

to  minimize  the  total  error,  one  has 
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(11) 


'  V, 

r  _ 

"  V, 

2;^  (A*,)2 

^ha(AvaAva) 

^u-t(AxtAzA) 

9x 

Z  wAv*  (<pk  -<p) 

4=1 

4=! 

4  =  1 

4=1 

AV 

v* 

^u'a(AyaA>'a) 

£wt(AytAzt) 

<Py 

~ 

2>,a,4(*-*> 

4=1 

4*1 

4  =t 

4=1 

N* 

*/> 

** 

V. 

^ha(AxaAza) 

Zwt(AzA)2 

<P_ 

ZnAz,^  -$o 

_  4  =  1 

4=1 

4=1 

4=1 

The  solution  to  the  liner  system  in  (1 1 )  gives  the  gradient  ealeulated  in  the  least-squares  sense. 
Convection  Schemes 

In  finite  volume  methods,  the  eonveetion  term  ean  he  calculated  using  the  Gauss  theorem 

|V(>VWQ=  \n<y<p)dY  (,2) 

n,  rr 

The  surface  integration  in  (12)  can  be  approximately  calculated  as 

\n{V<p)dV*Yjy-~S)f<pf  (13) 

r,  / 

The  convection  scheme  determines  how  the  face-centcr  value  <pt  is  calculated.  The  most  widely 

referenced  boundedness  criterion  using  the  normalized  variable  approach  is  Gaskell  and  Lau’s 
Conveetion-Boundedness  Criterion  (CBC)  (Gaskell  &  Lau4). 

Convection  Boundedness  Criterion  (CBC).  The  concept  of  CBC  is  easily  illustrated  in 
one-dimension  as  shown  in  Figure  4,  where  D  is  the  donor  cell,  U  is  the  upstream  eell,  and  A  is 
the  acceptor  cell. 


Figure  4.  Illustration  of  nodes  and  faee  in  onc-dimension 
Defining  a  normalized  variable 


<Pa~<Pu 

wc  thus  have 
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A,  =  *»-«'  ? 

<P.  -  Vu  and  S’, "«  ^ 

Based  on  the  normalized  variables,  CBC  states  the  following  local  boundedness  criterion 
<Pn  -  $7  -  1  for  e  [0, 1] 

^  04) 
Vf-Vo  for  $d€[ 0,1] 

The  graphical  representation  of  CBC  is  often  shown  in  the  normalized  variable  diagram  (NVD) 
(Leonard')  of  Figure  5. 


Figure  5.  Normalized  variable  diagram  (NVD)  with  CBC  (the  shaded  region) 


In  general,  the  normalized  face-center  value  can  be  written  as  a  function  of  the 
normalized  value  of  the  donor  cell 


9/  =f(9n) 


(15) 


Linear  Schemes.  If  the  function  /  in  (15)  is  linear,  the  seheme  is  called  a  linear  scheme. 
Examples  of  linear  schemes  include: 

Central  differencing  (CD)  scheme 


9f  = 


(16) 


Upwind  differencing  (UD)  scheme 


9/  =  9n  (17) 

Downwind  differencing  (DD)  scheme 

9/  =  1  (18) 


9 


Quadratic  Upstream  Interpolation  for  Convective  Kinetics  (QUICK)  scheme 


$f  -  ^  +  g  n  9) 

Warming  &  Beam  Second-Order  Upwind  (SOUD)  seheme 
*Pf  =  4  +  g  (2$) 

The  NVD  of  these  sehemes  and  the  CBC  region  are  shown  in  Figure  6. 


Figure  6.  NVD  of  linear  schemes  with  CBC 


Nonlinear  Schemes \  The  order  barrier  (Godunov6)  for  linear  schemes  implies  that  the 
CBC  sehemes  with  accuracy  of  sccond-ordcr  or  above  must  be  nonlinear  sehemes,  i.e.  the 
function  /  in  (15)  must  be  nonlinear.  Hereafter,  wc  summarize  some  of  the  CBC  schemes 
(limited  sehemes  in  NavyFOAM). 

van  Leer  Scheme  (van  Leer  ) 


j(<p)= 


2<p-(p\  ^e[0,l] 

<p,  ^£[0,1] 


(21) 


The  NVD  of  the  seheme  and  the  CBC  region  arc  shown  in  Figure  7. 
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Figure  7.  NVD  of  vanLeer  scheme  with  CBC 


Gamma  Scheme  (Jasak  et  al  ) 

~2 Yv  “(1  +  2^)^’ 


/(*)  = 


£e[0,/?] 

rc[o,i] 


with  —</?<- 
10  2 


The  NVD  of  the  scheme  and  the  CBC  region  are  shown  in  Figure  8. 


Figure  8.  NVD  of  Gamma  scheme  w  ith  CBC 


(22) 


IIYPER-C  Scheme  (Leonard) 


/{$)  = 


min(l,  -^(p),  ^e[0, 1] 
<p,  <P*  [0, 1] 


(23) 


The  HYPER-C  scheme  requires  that  the  local  Courant  number  C  f  <  1 .  The  NVD  of  the  scheme 
and  the  CBC  region  are  shown  in  Figure  9. 
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Figure  9.  NVD  of  HYPER-C  scheme  with  CBC 


Ultimate-Quickest  (UO)  Scheme  (Leonard  ) 

/(£)  = 


.  {%CfvH\-C(){b<p+l)  r 

min  {  g  ,  y  j 


<P , 


HYPER G  [0,  1] 

[0,1] 


(24) 


The  scheme  requires  that  the  loeal  Courant  number  C f  <  1 .  The  NVD  of  the  seheme  and  the 
CBC  region  are  shown  in  Figure  10. 


tCfj»(\-Cf)(6^3) 
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Compressive  Interface  Capturing  Scheme  for  Arbitrary  Meshes  (C1CSAM)  (Ubbink& 

Issa"l1 


f  )  —  Y j  f  \ HYPER  -C  )  "^  0  Y /)./" [IQ(<P  ) 


(25) 


1  +  cos  20, 


1 


0f  is  the  angle  between  the  normal  unit  veetor  of  the  front 


where  y  f  =  min 

V  ^  / 

(interface  between  two  phases)  and  the  vector  pointing  from  D  to  A,  see  Figure  1 1. 
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The  NVD  of  the  scheme  and  the  CBC  region  are  shown  in  Figure  12. 


i 

UQ 


Figure  12.  NVD  of  C1CSAM  scheme  with  CBC 
High  Resolution  Interface  Capturim;  Scheme  (HRICl  (Muzaferiia  &  Pencil  Let 


f\(p)  = 


2<p , 

L 

P, 


^6  [0,1] 
^e[i,l] 
pe  [0, 1] 


and 


Mp)  =  y/A(p)+0-y/)p 


with  yf  =  J\  cos  0f  |  , 


The  HRIC  scheme  can  be  written  as 


f(P)  = ' 


j f2(P )> 

°-^UP)4 


oi-cf 

04 


p, 


Cf  <  0.3 
0.3  <C/  <0.7 
Cf  >  0.7 


(26) 


The  NVD  of  the  scheme  and  the  CBC  region  are  shown  in  Figure  13. 
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Modified  HR1C  (MHEIQ  L Park  ct  ai /:)  Let 


2  <p. 


/(£)  =  U 


<p> 


^e[0,  j] 
^  e  [f ,  1] 
£e[0,l] 


and 


f2(<P)  = 


min(^,  ^e[0, 1] 

<P,  (P£[0, 1] 


MHR1C  can  be  written  as 

/(&)  =  r/f\(v)+Q-r/)Mv)  (27) 

with  yf  =  yj  |  cos  Qf  \  .  The  NVD  of  the  scheme  and  the  CBC  region  are  shown  in  Figure  1 4. 


Inter-Gamma  Scheme  (interGamma)  (Jasak  &  Weller 


-  2(p-  +  3 <p,  <p  e  [0, 


/(£)  =  <,! 


V, 


pe[0,l] 


(28) 
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The  NVD  of  the  scheme  and  the  CBC  region  are  shown  in  Figure  1 5. 


Figure  15.  NVD  of  inter-Gamma  scheme  with  CBC 


Modified  inter-Gamma  Scheme  (xnterGammaM)  Let 


-2(pl  +  3  <p, 

1 

<P. 


<P£{  OJ] 


be  the  original  inter-Gamma  scheme,  which  can  be  modified  as  follows 

f(v)  = 


Cf  <  0.3 


M<p), 

+  0.3 <Cf  < 0.7 


<P, 


Cf  >  0.7 


The  NVD  of  the  scheme  and  the  CBC  region  are  shown  in  Figure  16. 


Figure  16. 


NVD  of  interGammaM  scheme  with  CBC 


(29) 
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Modified  inter-Gamma  Scheme  (interGammaMD)  Let 


M9)  = 


-2<p2+3p,  £e[0,j] 

9  &[\i  i] 


1 

9, 


be  the  onginal  inter-Gamma  scheme,  and 

Uv)-r,m*U-r,->¥  Wiihr/-Vi^ 

The  interGammaMD  scheme  can  be  written  as 
f 2(9)1 


f(9)  = 


Cf  <  0.3 


^  fi (9)  +  (l  -  ^5%  0.3  <  Cf  <  0.7 


9i 


Cf  >  0.7 


(30) 


(31) 


Interpolation  Schemes 

Cell-based  Surface  Interpolation  Schemes , 

In  cell-centcr-based  finite  volume  methods,  it  is  often  required  to  calculate  the  faee- 
centcr  value.  The  interpolation  schemes  needed  for  non-convection  terms  will  be  introduced  in 
this  section.  The  situation  is  illustrated  in  Figure  17.  The  owner  and  neighbor  cells  of  the  face 
may  or  may  not  be  located  within  a  mesh  block  of  a  single  processor. 

\ 

Owner  processor  \  Neighbor  processor 

\  Processor  interface 


Figure  17.  Cell-based  interpolation  of  face-center  value 


Linear  Interpolation  Surface  Interpolation  Scheme .  The  linear  interpolation  calculates 
the  face-ccnter  value  as 

<Pf  =  A<pp  +  (1  -  X)q>N  (32) 

where  A  is  the  weighting  factor  calculated  as 
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The  linear  interpolation  scheme  may  produce  large  errors  due  to  non-orthogonality  and  skewness 
of  unstructured  meshes. 

reconCentral  Surface  Interpolation  Scheme .  The  face-centcr  value  can  he  reconstructed 
using  both  the  value  and  the  gradient  at  neighboring  cell  centers. 

<P/  =^k>  +  (V v)r  ■  ?pf  +  +  (V <P)N  '  rv]  (33) 


Upwind  Deferred  Correction  (UPDC)  Surface  Interpolation  Schemes.  To  improve 
stability,  the  face-centcr  value  can  sometimes  be  calculated  using  the  concept  of  deferred 
correction  (Khosla  &  Rubin14;  Hayase  et  al.15). 

9f-9f  +{<Pf-<Pf  )  (34) 

where  <pv,OK  is  the  value  calculated  using  a  first-order  upwind  scheme,  and  (px)  is  the  value 

calculated  using  higher-order  interpolation  schemes.  The  superscript  “old”  represents  the 
previous  time  or  iteration  step.  The  candidates  for  the  higher-order  scheme  may  include  the 
reconCentral  scheme  and  some  of  the  higher-order  limited  schemes  in  the  next  section. 

The  first-order  upwind  scheme  can  be  written  as 

<Pf°l  =sign+(F/)pP  +  [l-sign+(F/)]$?v  (35) 


where  Ff  =  nf  •  V f  is  the  volume  flux  through  the  face,  and 


signYF,;  = 


Ff  >0 
Ff  <0 


The  higher-order  normalized  variable  based  limited  schemes  can  be  written  as 


f  —  typ  ^ typ) 


where  'F(r)  is  the  limiter  function.  Substituting  (35)  and  (36)  into  (34)  yields 


(36) 


<Pf  = 


<Pp  +  [^(r)(<Ps  ~  <Pp)]°Ui 

i 


Ff  >0 

Ff  <0 


(37) 
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Node-Based  Surface  Interpolation  Schemes . 

We  consider  a  face  of  any  polyhedral  cell 


Figure  18.  Node-based  interpolation  of  faee-eenter  value 

Let  Nv  be  the  total  number  of  vertices  of  faces,  the  faee-eenter  value  is  calculated  as  the  average 
of  the  nodal  values 

1  y 

TZX  (38) 


The  nodal  values  are  calculated  using  volume-to-point  interpolation  described  in  the  following 
section. 

Volume-to-Point  Interpolation  Schemes 

It  is  easier  to  illustrate  in  two-dimensions  and  the  extension  to  three-dimension  is  rather 
straightforward. 


Figure  19.  Stencil  in  cell-based  gradient  calculation  in  two-dimensions 


In  volume-to-point  interpolation,  the  nodal  value  cpn  is  calculated  as  a  weighted 
averaging  of  surrounding  eell-eenter  values 

Nm  Nn 

<Pn='E(<PcJWcJ/'ZWcj)  (39) 

/=!  j= I 


where  Nn  is  the  total  number  of  neighboring  cells  of  node  //,  wci  is  the  weighting  factor. 

V o lu m e -to- Point  Interpolation  Based  on  PLWA .  In  P seudo- Laplacian- Weight ed- 
Averaging  (PLWA)  (Holmes  &  Connell  6;  Frink  ;  Kim  et  al.  s),  the  weighting  factors  in 
Equation  (39)  are  calculated  by  solving  the  following  optimization  problem. 
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Giving  the  constraint 

L(Xr,)  =  ^WcAXc.,-Xn)=() 

1=1 

Uy„)  =  YdwCJ(yeJ-yn)  =  0  (40) 

1  =  1 

»=l 

we  need  to  find  the  weighting  factors  wc  i  that  minimize  the  cost  function 

C  =  'Z(rc,Awc.)2  (41) 


r- 1 


where  rci  =  (xc yci ,  zci )  is  the  position  vector  of  the  cell  center,  rn  =  (xn ,  yn  ,zn)  is  the  position 

vector  of  the  node  n,  and  rCJ  =|  rcj  -  Fn  \=  ^(xci  -  .r„  )2  +  (yc,  -  y„  )2  +  (=CJ  ~=„Y  .  The  Aw, ,  is 
related  to  the  weighting  factor  by  wcJ  =  1  +  Aw. . . 

Using  the  method  of  Lagrange  multipliers,  the  Lagrange  function  is  defined  as 

A(  wr , ,  we,: we  v> ;  Ax ,  Ay,  A  )  =  C- 2[AxL(.xn )  +  AvL(yn )  +  A.L(zH )]  (42) 

The  optimal  solution  is  found  by  solving  the  following  equation 

V,Vl  . „■  Wc.2 . WcM.  *  '  'O  =  0  (43) 

which  can  be  written  in  matrix  form 

r  0 1 

(44) 

[Zirj  [v\  [a  L"  tv 

where 


[diag(/-2)] 

-[Ar] 

Aw 

0  ‘ 

.  [Ar]T 

[0]  . 

„  k 

-R 

[diagfr2 )]  = 


rc.l 


c.2 
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[Ar]  = 


*,..-■**  yc.i-y»  2c.\~zn 

Xc.2~Xn  yc.2~y*  Zc.2  ~  2n 
Xc.N.~\  Xt.Af,  ~yn  Zc.N.  ~  =n 


Aw  = 


Am’ 


c.l 


Aw. 


c.2 


Avv 


c.M, 


X  = 


and 


R  = 

R, 

- 

-jo 

R 

Solving  the  linear  system  of  equations  we  obtain 

*=- nr'R  and 

(45) 

Aw  =  [diag(r:)]~'[Ar]X 

where 

[I]  =  [Ar]  T[diag(r2)]^l[Ar] 


Volume-to-Point  Interpolation  Based  on  IDWA.  In  Inverse-Distance- Weighted- 
Averaging  (IDWA),  the  weighting  factors  in  Equation  (39)  are  caleulated  based  on  the  distance 
between  the  node  and  each  neighboring  cell  center. 


with 


r  —\r  —  r 

cj  I  cj  n 


1=  \l(xc.,-x„)2  +(yc,~y„)2  +(zc 
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Diffusion  Schemes 

In  finite  volume  methods,  the  diffusion  term  can  be  calculated  using  the  Gauss  theorem 


j V  {yS/(p)dQ.  =  \h  (yV(p)dr 

a  r. 


(47) 


The  surface  integration  in  (47)  can  be  approximately  calculated  as 


\n(yV<p)dY  «  •  V(p)fSf  =X^/5/) 


d(p 
\dn  J 


(48) 


Figure  20.  Calculating  face-normal  gradient 
The  face-normal  gradient  can  be  calculated  as 

f 


d<p 
dn  j. 


=  T 


Ar, 


/ 

\ - v - * 

snGrad  correchon 


(49) 


where  Cf  =  ///  -  A rrN  /  |  A rPN  |  is  the  non-orthogonal  correction  vector,  and 

¥*>,  =  W  <pP  +  (1  -  A)V  <pN 
A  is  the  inverse-distance  weighting  factor. 

Solution  Algorithms 

Pressure-Velocity  Coupling 

Solving  the  Pressure  Equation .  The  momentum  equation  can  be  written  in  discretized 
form  as  a  system  of  linear  equations 

AF’+V//  =7  (50) 

where  V*  is  the  velocity  vector. 
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p  is  the  guessed  pressure. 


Vp  = 


n 


A  =  D-B  =  (D„,+DJ-B  is  the  matrix  of  the  linear  equations, 
D  =  (D.„  +  )  =  diag(A) 


Din  =  part  of  diag(A)  contributed  from  internal  faces 
D^.  =  part  of  diag(  A)  contributed  from  boundary  faces 
—  B  =  A  —  D  =  off-diagonal  part  of  A, 

f  =  fin  4  fhc  ~  the  source  vector  absorbing  any  explicit  term  and  source  term, 
fin  =part  of  source  vector  contributed  from  internal  faces,  and 

ft*  =Part  of  source  vector  contributed  from  boundary  faces. 

Equation  (50)  can  be  written  as 


(D-B  )V'+Vp'=f 


(51) 


Giving  guessed  pressure,/?*,  Equation  (51)  is  solved  for  velocity,  V* .  This  is  the 
predictor  step  of  the  SIMPLE  or  PISO  method.  Because  the  velocity  obtained  from  the  predictor 
step  doesn’t  satisfy  the  continuity  equation,  both  the  pressure/?  and  velocity  V*  need  to  be 
corrected. 

In  order  to  derive  the  discretized  pressure  equation,  the  following  algebraic  manipulation 
was  applied  to  the  momentum  equation.  The  matrix  A  of  discretized  momentum  equations  can 
be  written  as 


(52) 


with 


Dcrnpfav  =  component  -average  part  of  diag( A) contributed  from  boundary  faces  . 


Let 


and 
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thus 


A  =  (D„  +  Djr" )  -  (B  +  DT"*  -  D* )  =  Dd  -  Bn 
Now,  Equation  (50)  can  be  written  as 

(D0-B  D)V'  +  Vp=J  (53) 

Let  v  and  p"  be  corrected  velocity  and  pressure,  respectively.  They  should  satisfy  the 
discretized  momentum  equation  (53),  i.e. 

(D0-B0)r  +  v=7  <54) 

which  can  be  approximated  as 

D0K*’-B  DV'+Vp'=f  (55) 

Because  Dd  in  (55)  is  a  diagonal  matrix,  it  is  trivial  to  solve  (55)  to  obtain 

v "  =  D-' (B0  P‘  +  7)  -  D-' V*  =  P  ~  vfip'  (56) 

where  V  =  D~'  (B/}  V*  +  f)  is  the  pseudo-velocity  vector. 

Substituting  (56)  into  the  continuity  equation,  V  •  V  =0,  we  obtain  the  following  discretized 
Poisson  equation  for  pressure 

V-(D-'V//*)  =  V-P  (57) 

It  is  shown  later  that  Equation  (57)  incorporates  the  idea  of  the  Rhic-Chow  momentum 
interpolation  scheme. 

In  the  PI  SO  method,  consecutive  corrector  steps  may  be  used  to  correct  pressure  and 
velocity,  the  momentum  equation  that  is  satisfied  at  the  (A+l)-th  steps  is 

D„  K*+l  -  B„  Vk  +  Vp*+I  =  f  (58) 

thus 

Vk"  =Do'(B0E*  +  /)-I)‘lV//t|  -  V  -  DpV/?<+l  (59) 

where  V  =  DJ,1  (B„  Vk  +  f) .  Note  that  k  =  0  represents  the  predictor  step,  i.e.  V°  =  V' . 
Substituting  (59)  into  the  continuity  equation  V  •  Vktl  =  0  yields 

V  •  (DpVpul)  =  V  P  (60) 
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User’s  Guide 


This  section  provides  information  for  users  to  help  with  running  any  of  the  updates 
developed  specifically  as  a  part  of  NavyFOAM  V1.0.  In  addition,  there  is  a  supplement  to  the 
original  OpenFOAM  User’s  Guide  contained  in  Appendix  A  that  could  be  useful  to  readers.  For 
more  detailed  information  on  the  OpenFOAM  code  and  settings  consult  the  OpenFOAM  User’s 
Guide:  http://foam.sourceforge.net/doc/Guides-a4/UserGuide.pdf.  In  addition,  to  aid  users 
several  utilities  have  been  developed  to  aid  in  post-processing  NavyFOAM  results,  which  are 
described  in  Appendix  B.  Finally,  tutorials  have  been  developed  that  aid  a  user  in  running  the 
codes  which  include:  a  lid  driven  cavity  problem,  a  fully  submerged  axisymmetric  body  and  the 
Wigley  hull  with  free  surface  in  Appendices  C,  D  and  E,  respectively.  It  is  recommended  that 
new  user’s  go  through  the  tutorials  to  learn  the  preferred  settings  to  use  with  the  solvers. 

Gradient  Schemes 

Navy  Least  Squares  Gradient  Schemes 

The  gradient  scheme  is  specified  by  a  sub-dictionary  entry  gradSchemes  in  the  system 
finite  volume  dictionary  file  fvSchemes.  Users  may  refer  to  the  OpenFOAM  User  Guide  Ver.  1.5 
Section  4.4.3  on  page  U-l  10  for  more  details  about  this  sub-dictionary.  Figure  21  illustrates  an 
abbreviated  example.  The  default  gradient  scheme  is  the  one  using  Gauss  theorem  with  the  face- 
center  value  calculated  by  linear  interpolation.  The  scheme  used  in  calculating  the  gradient  of 
pressure  is  NavyLeastSquares  -  where  the  gradient  is  calculated  in  a  least  squares  sense,  more 
details  regarding  the  NavyLeastSquares  gradient  scheme  can  be  found  in  the  Technical 
Description  section 

Dictionary  file: 

$  (CASE_DIR)  / system/ fvSchemes 

Sub-dictionary: 

gradS  cheme  s 

{ 

//-  default  scheme 
default  Gauss  linear; 

//-  gradient  scheme  for  pressure  p 
grad (p)  NavyLeastSquares ; 

} 

Figure  21.  Example  of  gradSchemes  sub-dictionary 
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Interpolation  Schemes 

Cell-based  Surface  Interpolation  Schemes 

reconCentral  Surface  Interpolation  Scheme  The  reconCentral  surface  interpolation 
scheme  calculates  the  face-center  value  of  <f>  using  both  the  value  and  the  gradient  of  <f>  at  the 
center  of  the  owner  and  neighbor  cells.  The  following  example  in  Figure  22  shows  how  to 
specify  the  interpolation  schemes  in  the  sub-dictionary  entry  interpolationSchemes  of  the 
system  finite  volume  dictionary  file  fvSchemes.  Users  may  refer  to  OpenFOAM  User  Guide  Ver. 
1.5  Section  4.4.1  on  page  U-108  for  more  details  about  this  sub-dictionary'.  The  default  scheme  is 
linear  interpolation,  where  the  face-center  value  is  calculated  as  a  weighted  average  of  the  value 
at  the  center  of  the  owner  and  neighbor  cells.  The  weighting  factors  are  based  on  inverse 
distance.  The  reconCentral  scheme  is  used  to  calculate  the  face-center  value  of  u. 


Dictionary  file: 

$  (CASE_DIR)  /system/ fvSchemes 

Sub-dictionary: 

interpolationSchemes 

{ 

//-  default  scheme  is  linear  interpolation 
default  linear; 

//-  surface  interpolation  of  U 
interpolate (U)  reconCentral ; 

} 

Figure  22.  Example  of  interpolationSchemes  sub-dictionary 

reconCentralDC  Surface  Interpolation  Scheme  The  reconCentral  is  almost  a  second- 
order  interpolation  scheme.  Besides  calculating  the  surface  interpolation  at  the  facc-center,  it 
may  also  be  used  to  improve  the  accuracy  in  discretizing  the  convection  term  in  the  momentum 
equation  using  deferred  correction.  The  following  example  in  Figure  23  shows  the  sub-dictionary 
divSchemes  in  the  finite  volume  system  file  fvSchemes .  Users  may  refer  to  OpenFOAM  User 
Guide  Ver.  1.5  Section  4.4.5  on  page  U-l  1 1  for  more  details  about  this  sub-dictionary.  The 
convection  term  is  integrated  using  the  Gauss  theorem  and  the  face-ccnter  velocity  U  in  the 
convection  term  takes  the  Upwind  Deferred  Correction  (UPDC)  form  using  the  reconCentral 
scheme  to  calculate  the  higher-order  correction.  More  details  regarding  UPDC  can  be  found  in 
the  Technical  Description  Section. 
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Dictionary  file: 

$  (CASE__DIR)  /system/fv Schemes 

Sub-dictionary; 

div Schemes 

{ 

//-  convection  term  in  momentum  equation 
div (phi,  U)  Gauss  reconCentralDC; 

> 

Figure  23.  Example  of  divSchemes  sub-dictionary 
Node-based  Surface  Interpolation  Schemes 

reconPLWA  Surface  Interpolation  Scheme  The  reconPLWA  surface  interpolation 
scheme  calculates  the  face-center  value  as  an  average  of  the  nodal  values  on  the  face.  The  nodal 
value  is  calculated  from  cell-center  values  using  a  volume-to-point  interpolation  scheme  based 
on  Pseudo-Laplacian-Weighted-Averaging  (PLWA).  More  details  of  PLWA  can  be  found  in  the 
Technical  Description.  The  following  example  in  Figure  24  shows  the  sub-dictionary  entry 
interpolationSchemes  in  the  system  finite  volume  dictionary  file  fvSchemes.  Users  may 
refer  to  OpenFOAM  User  Guide  Ver.  1.5  Section  4.4.1  on  page  U-108  for  more  details  about  this 
sub-dictionary.  The  default  interpolation  scheme  is  linear.  The  next  line  of  the  dictionary  shows 
that  the  reconPLWA  scheme  is  used  to  calculate  the  face-ccntcr  value  of  U 

Dictionary  file; 

$  (CASE_DIR)  /system/fvSchemes 

Sub-dictionary; 

interpolationSchemes 

{ 

//-  default  scheme  is  linear  interpolation 
default  linear; 

//-  surface  interpolation  of  U 
interpolate (U)  reconPLWA; 

> 

Figure  24.  Example  of  interpolationSchemes  sub-dictionary' 

reconlDWA  Surface  Interpolation  Scheme  The  reconIDWA  surface  interpolation 
scheme  calculates  the  faee-eenter  value  as  an  average  of  the  nodal  value  on  the  face.  The  nodal 
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value  is  calculated  from  cell-center  values  using  a  volume-to-point  interpolation  scheme  based 
on  Inverse-Distance-Weighted-Averaging  (IDWA).  More  details  of  IDWA  can  be  found  in  the 
Technical  Description.  The  following  example  in  Figure  25  shows  the  sub-dietionary  entry 
interpolationSchemes  in  the  system  finite  volume  dictionary  file  fvSchemes.  Users  may 
refer  to  OpenFOAM  User  Guide  Ver.  1.5  Section  4.4.1  on  page  U-108  for  more  details  about  this 
sub-dictionary.  The  default  interpolation  scheme  is  linear.  The  reconiDWA  scheme  is  used  to 
calculate  the  face-center  value  of  u. 


Dictionary  file: 

$  (CASE_DIR)  /system/fvSchemes 

Sub-dictionary: 

interpolationSchemes 

{ 

//-  default  scheme  is  linear  interpolation 
default  linear ; 

//-  surface  interpolation  of  U 
interpolate (U)  reconiDWA; 

) 

Figure  25.  Example  of  interpolationSchemes  sub-dietionary 

Convection  Schemes 

Convection  Boundedness  Criterion  (CBC)  Schemes 

The  following  NVD  ( normalized  variable  diagram)  based  CBC  schemes  are  particularly 
useful  in  capturing  interlaces  between  two  fluids,  e.g.  air  and  water,  in  the  two-phase  flow  solver 
of  NavyFOAM.  More  details  of  CBC  schemes  can  be  found  in  the  Technical  Description. 

Compressive  Interface  Capturing  Scheme  for  Arbitrary  Meshes  (C ICS  AM)  The 
following  example  in  Figure  26  shows  the  sub-dictionary  entry  divSchemes  of  the  system  finite 
volume  dictionary  file  fvSchemes  that  specifies  the  divergence  schemes  for  the  convection  term 
in  the  transport  equation  of  volume  fraction  (gamma).  In  this  example  the  convection  term  is 
integrated  using  the  Gauss  theorem  and  the  cicsam  seheme  is  used  to  calculate  the  face-center 
value  of  gamma.  Users  may  refer  to  OpenFOAM  User  Guide  Ver.  1.5  Section  4.4.5  on  page  U- 
111  for  more  details  about  this  sub-dictionary. 
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High  Resolution  Interface  Capturing  Scheme  (HRIC)  The  following  example  in 
Figure  27  shows  the  sub-dictionary  entry  divSchemes  of  the  system  Finite  volume  dictionary 
file  fvSchemes  that  specifies  the  divergence  schemes  for  the  convection  term  in  the  transport 
equation  of  volume  fraction  (gamma).  In  this  example,  the  convection  term  is  integrated  using 
the  Gauss  theorem  and  the  hric  scheme  is  used  to  calculate  the  face-center  value  of  gamma. 
Users  may  refer  to  OpenFOAM  User  Guide  Ver.  1.5  Section  4.4.5  on  page  U-l  1 1  for  more 
details  about  this  sub-dictionary. 


Modified  HRIC  (MHRIC)  The  following  example  in  Figure  28  shows  the  sub-dictionary 
entry  divSchemes  of  the  system  finite  volume  dictionary  file  fvSchemes  that  specifies  the 
divergence  schemes  for  the  convection  term  in  the  transport  equation  of  volume  fraction 
(gamma).  In  this  example,  the  convection  term  is  integrated  using  the  Gauss  theorem  and  the 
MHRIC  scheme  is  used  to  calculate  the  face-center  value  of  gamma.  Users  may  refer  to 
OpenFOAM  User  Guide  Ver .  1.5  Section  4.4.5  on  page  U-l  1 1  for  more  details  about  this  sub¬ 
dictionary. 
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Inter-Gamma  Scheme  (interGamma)  The  following  example  in  Figure  29  shows  the 
sub-dictionary  entry  divSchemes  of  the  system  finite  volume  dictionary  file  fvSchemes  that 
specifies  the  divergence  schemes  for  the  convection  term  in  the  transport  equation  of  volume 
fraction  (gamma).  In  this  example,  the  convection  term  is  integrated  using  the  Gauss  theorem 
and  the  interGamma  scheme  is  used  to  calculate  the  face-center  value  of  gamma.  Users  may 
refer  to  OpenFOAM  User  Guide  Ver.  1.5  Section  4.4.5  on  page  U- 111  for  more  details  about  this 
sub-dietionary. 


Dictionary  file: 

$  (CASE_DIR)  / system/fvSchemes 

Sub-dictionary: 

divSchemes 

{ 

//-  convection  term  in  gamma  equation 
div (phi , gamma)  Gauss  interGamma; 

} 

Figure  29.  Example  of  divSchemes  sub-dietionary 

Modified  inter-Gamma  Scheme  (interGammaM)  The  following  example  in  Figure  30 
shows  the  sub-dictionary  entry  divSchemes  of  the  system  finite  volume  dictionary  file 
fvSchemes  that  specifics  the  divergence  schemes  for  the  eonveetion  term  in  the  transport 
equation  of  volume  fraction  (gamma).  In  this  example,  the  convection  term  is  integrated  using 
the  Gauss  theorem  and  the  interGammaM  scheme  is  used  to  calculate  the  faee-eentcr  value  of 
gamma.  Users  may  refer  to  OpenFOAM  User  Guide  Ver.  1.5  Section  4.4.5  on  page  U-l  1 1  for 
more  details  about  this  sub-dictionary. 


i 
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Dictionary  file: 

$  (CASE_DIR)  /system/fvSch ernes 

Sub-dictionary: 

divSch ernes 

{ 

//-  convection  term  in  gamma  equation 
div (phi , gamma)  Gauss  interGammaM; 

} 

Figure  30.  Example  of  divSchemes  sub-dictionary 

Modified  Inter-Gamma  Scheme  (interGammaMD)  The  following  example  in  Figure 
31  shows  the  sub-dictionary  entry  divSchemes  of  the  system  Finite  volume  dictionary  file 
fvSchemes  that  specifies  the  divergence  schemes  for  the  convection  term  in  the  transport 
equation  of  volume  fraction  (gamma).  In  this  example,  the  convection  term  is  integrated  using 
the  Gauss  theorem  and  the  interGammaMD  scheme  is  used  to  calculate  the  face-ccnter  value  of 
gamma.  Users  may  refer  to  OpenFOAM  User  Guide  Ven  i.5  Section  4.4.5  on  page  U-l  1 1  for 
more  details  about  this  sub-dictionary. 

Dictionary  file: 

$  (CASE_DIR)  /system/fvSchemes 

Sub-dictionary: 

divSchemes 

{ 

//-  convection  term  in  gamma  equation 
div (phi , gamma)  Gauss  interGammaMD ; 

} 

Figure  31.  Example  of  divSchemes  sub-dictionary 
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Example  Results 

Results  are  demonstrated  for  a  fully  submerged  axisymmetrie  body,  Body-1,  the 
K.VLCC2  tanker  without  a  free  surfaee  using  the  double  body  approximation,  DTMB  Model 
5415  with  both  fixed  sinkage  and  trim  as  well  as  dynamie  sinkage  and  trim  and  the  Joint  High 
Speed  Sealift  (JHSS)  eoneept  surface  ship  with  and  without  wateqet  propulsion. 

Body-1 

This  seetion  involves  using  NavyFOAM's  steady,  incompressible  Reynolds  Averaged 
Navier-Stokes  (RANS)  solver  for  a  3-D  body-of-revolution  referred  to  as  Body-1.  The  RANS 
equations  are  solved  using  NavyFOAM’s  k-omega  SST  turbulence  model.  Only  half  the  body  is 
solved,  as  symmetry  is  assumed,  and  the  domain  is  non-dimensionalized  by  length.  The 
Reynolds  number  (Re)  based  upon  the  body  length  is  6.6  million.  The  boundary  conditions  used 
for  these  computations  are: 

•  Defined  fixed  turbulent  quantities  (k,  omega,  nuTilda)  and  velocity  (U)  at  inlet 

•  Defined  pressure  (p)  at  outlet 

•  Zero  gradient  for  all  quantities  at  farfield  boundaries 

•  nuTilda  =  0,  k  ~  0,  and  omega  set  to  zero  gradient  at  the  walls 


A  side  view  of  the  ONR  Body  1  geometry  can  be  seen  below  in  Figure  32. 


Figure  32.  Body- 1  geometry 


Computations  were  done  on  unstructured  meshes  that  eontain  tetrahedral  (in  the  flow 
field)  and  prism  (in  the  boundary  layer)  elements.  Grids  for  these  computations  typically 
contained  approximately  2  million  eells  total.  The  left  of  Figure  33  shows  an  unstructured 
surface  mesh  on  the  body  and  the  symmetry  plane,  and  the  right  side  of  Figure  33  shows  a 
surface  mesh  w  ith  a  span-wise  cross  sectional  cut  of  the  volume  mesh  at  the  midbody. 
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Figure  33.  Surface  mesh  on  the  body  and  symmetry  plane  (left)  and  surface  mesh  with  volume 

mesh  cross  cut  (right) 


The  computational  domain  is  split  into  many  domains  to  allow  the  computations  to  be  run 
in  parallel.  The  domain  is  split  using  the  METIS  domain  decomposition  method.  Typical  steady 
state  run  times  for  this  geometry  are  3-5  hours  depending  on  the  number  of  domain  partitions. 
Steady  state  convergence  is  assumed  when  the  forces  (pressure  and  viscous)  on  the  body  change 
by  a  negligible  amount  from  one  iteration  to  the  next.  Figure  34  shows  results  from  NavyFOAM 
computations. 


Figure  34.  Axial  velocity  (Ux)  contours  on  the  symmetry  plane  and  pressure  (Press)  contours  on 

the  hull 

One  can  notice  that  the  stagnation  point  is  qualitatively  predicted  well  at  the  nose  of  the 
body.  The  velocity  slows  to  zero  and  the  pressure  on  the  body  is  a  maximum  at  the  nose.  While 
the  drop  in  pressure  associated  with  the  acceleration  of  the  fluid  around  the  shoulder  of  the  bow 
is  also  predicted  correctly.  Computations  also  show  the  flow  remaining  attached  along  the  body 
leaving  a  wake  after  the  stem. 

Figure  35  shows  some  of  the  quantitative  results  from  the  NavyFOAM  computations 
compared  to  experimental  measurements. 
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Figure  35.  Skin  friction  coefficient  (Q)  and  Pressure  eoeffieient  ( Cp )  plotted  along  the  length  of 

the  body 

In  Figure  35  there  is  a  slight  disagreement  in  predieted  and  measured  Qat  the  tail  of  the 
body,  but  the  trends  are  matched  well.  The  NavyFOAM  Cp  predictions  match  experimental 
measurements  well  along  the  body. 

Figure  36  below  shows  comparisons  of  computed  and  measured  boundary  layer  profiles 
at  various  locations  along  the  body. 


Figure  36.  Axial  velocity  boundary  layer  plots  at  x/L  =  0.755  (left),  x/L  =  0.846  (middle),  x/L  = 

0.934  (right) 

The  axial  velocities  are  non-dimensionalized  by  the  free  stream  velocity  and  boundary 
layer  lengths  are  non-dimensionalized  by  the  radius.  NavyFOAM  computed  boundary  layer 
velocity  profiles  match  experimental  measurements  very  well  at  various  locations  along  the 
body. 

In  conclusion,  this  effort  shows  that  NavyFOAM  has  the  capability'  to  successfully 
predict  quantitative  and  qualitative  characteristics  for  a  body-of-revolution,  as  demonstrated  on 
the  Body-1.  The  RANS  computations  were  carried  out  successfully  on  multiple  processors  on 
unstructured  meshes.  Results  for  non-bodv-of-revolution  hull  forms  can  be  found  in  Delaney  et 
al.19. 
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KVLCC2 


The  KVLCC2  is  a  tanker  used  for  eode  validation  by  many  organizations  and  as  a  part  of 
several  international  workshops.  The  stem  of  the  KVLCC2  with  predicted  streamlines  is  show  n 
in  Figure  37,  flow  is  from  left  to  right. 


Figure  37.  Stem  flow  of  the  K.VLCC2 

We  ran  sRansFOAM  for  this  double-body  flow  case  in  which  the  free  surface  is  replaced 
by  a  symmetry  plane.  We  deliberately  ehose  to  use  two  unstructured  meshes  to  evaluate  spatial 
aeeuraey  of  the  predictions.  One  of  them,  shown  at  the  left  of  Figure  38,  is  a  hybrid 
unstructured  mesh  generated  using  SolidMesh  developed  at  the  SimCenter  of  the  Mississippi 
State  University.  The  mesh  consists  of  prisms  with  triangular  bases  near  the  hull  surfaec  and 
tetrahedral  in  the  rest  of  the  domain.  The  total  eell  eount  is  approximately  8  million.  The  other 
mesh,  shown  at  the  right  of  Figure  38,  is  a  hexahedron-dominant  unstructured  mesh  generated 
using  the  snappyHexMesh  utility  available  in  the  standard  OpenFOAM  paekage.  With  the  latter, 
we  put  in  three-levels  of  embedded  fine-mesh  bloeks  to  better  resolve  the  near-body  region.  Its 
total  number  of  elements  is  a  little  over  3  million.  The  near-wall  mesh  resolutions  for  both 
meshes  are  sueh  that  the  distance  front  the  wall  is  larger  than  40  wall  units  (y+  >  40)  over  a  large 
portion  of  the  hull  surface.  So,  the  wall  function  approach  was  employed  to  provide  the 
boundary  conditions  for  the  momentum  equations  and  the  turbulence  equations.  Wileox’s  k-O) 
model  (Wilcox1)  was  used  for  turbulenee  elosure  for  its  good  track  reeord  for  this  elass  of  flows 
(Kim  and  Rhee"')). 


Figure  38.  Grids  used  for  the  KVLCC2 

The  contours  of  mean  axial  velocity  ( U)  at  the  propeller  plane  (x/L  =  0.9823)  are  shown 
in  Figure  39  for  the  two  meshes.  As  ean  be  seen,  the  characteristic  shape  of  the  U-eontours 
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(“hook”-shaped  or  “rabbit’s”  ear-like)  is  closely  captured  by  both  meshes.  The  contours  of 
turbulent  kinetic  energy  at  the  same  plane  are  depicted  in  Figure  40.  The  region  of  high 
turbulent  kinetic  energy,  which  originates  from  the  upstream  boundary  layer  and  overlaps  with 
the  region  of  low  axial  velocity  depicted  in  Figure  39  is  reproduced  reasonably  well  by  the 
computations,  although  their  peak  values  are  somewhat  under  predicted.  The  overall  prediction 
accuracy  shown  here  with  the  unstructured  meshes  used  in  this  study  is  remarkably  good.  More 
details  on  these  calculations  can  be  found  in  Kim  et  al.  '. 


Figure  39.  Contour  of  axial  velocity  at  x/L  =  0.9825  predicted  on  the  two  meshes.  Top  hybrid 
(prism  +  tet)  unstructured  mesh;  Bottom  -  snappyHcxMesh 


Figure  40.  Contour  of  turbulent  kinetic  energy  at  x/L  =  0.9825  predicted  on  the  two  meshes.  Top 
-  hybrid  unstructured  (prism  +  tet);  Bottom  -  snappyHcxMesh 

DTMB  Model  5415 

Fixed  Sinkage  and  Trim 

These  computations  were  carried  out  using  ransFSFOAM  for  a  single  Froude  number  of  0.28 
using  three  different  isotropic  eddy-viscosity  based  turbulence  models  including  the  realizable  k- 
£■  (Shih  ct  al."),  SST  k-oj  (Mcnter*)  and  Wilcox’s  k-co  (Wilcox3)  models.  The  computational 
meshes  were  generated  using  GridPro  (www.undpro.com).  a  commercial  meshing  package  well 
known  for  high-quality  hexahedral  meshes.  Great  care  was  taken  to  ensure  that  such  salient 
features  as  the  hull-generated  waves,  the  boundary  layer  along  the  entire  hull,  and  the  near-wake 
are  adequately  resolved.  To  check  grid-dependency  of  the  solutions,  three  systematically  refined 
hexahedral  meshes  were  used  with  13  million  (fine),  6  million  (medium),  and  3  million  (coarse) 
elements.  One  of  the  grids  used  is  show  n  in  Figure  41. 
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Figure  41.  Gnd  for  DTMB  Model  5415 


Although  a  time-marehing,  transient  solution  algorithm  is  employed  in  ransFSFOAM, 
our  goal  was  to  obtain  steady-state  solutions.  An  iterative  implicit  solution  algorithm  employed 
in  ransFSFOAM  greatly  accelerated  solution  convergence  by  allowing  us  to  use  a  fairly  large 
time-step  size.  The  transient  computations  were  continued  for  sufficiently  long  periods  of  time 
until  not  only  the  global  quantities,  such  as  the  forces  and  moments  acting  on  the  hull,  but  also 
other  flow  features  like  wave  elevation,  hull  boundary  layer  and  wake,  all  reach  steady  states. 

Figure  42  illustrates  the  wave  pattern  predicted  with  the  SST  k-(0  model  result  on  the  6 
million  cell  mesh,  along  with  the  measured  one.  Figure  43  shows  the  wave  elevations  along  the 
three  longitudinal  cuts  (y/L  =  0.082,  0.172,  0.301)  obtained  using  the  SST  k-a)  model  on  all  three 
meshes.  First  of  all,  the  results  indicate  grid-eonvergenee  of  the  predictions,  all  of  which  are  in 
excellent  agreement  with  the  data.  The  results  obtained  using  three  different  turbulence  models 
on  the  medium  (6  million  cell)  mesh  are  shown  in  Figure  44  The  differences  among  the  three 
results  are  measurable  yet  small.  The  realizable  k-e  model  results  appreciably  deviate  from  the 
other  two  k-co  model  results. 
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Figure  42.  Contour  of  wave  elevation  for  DTMB  54 1 5  with  SST  k-a)  model  result  on  the  6 

million  cell  mesh 
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Figure  43.  Wave  elevations  along  three  longitudinal  cuts  obtained  using  SST  k-(0  model  on  three 
different  meshes:  top  -  y/L  =  0.082;  middle  -  y/L  =  0. 1 72;  bottom  -  y/L  =  0.30 1 


Figure  44.  Wave  elevations  along  three  longitudinal  euts  obtained  using  three  different 
turbulence  models  on  the  6  million  cell  mesh  :  top  -  y/L  =  0.082;  middle  -  y/L  =  0. 1 72;  bottom 

y/L  =  0.301 


Figure  45  shows  the  contours  of  axial  velocity  at  x/L  =  0.935  obtained  using  the  three 
turbulence  models  on  the  6  million  cell  mesh.  All  three  turbulence  models  eapture  the  gross 
feature  of  the  boundary  layer  -  the  distorted  velocity  contours  reflecting  thickening  of  the  stem 
boundary  layer  due  to  convergence  of  wall-limiting  streamlines,  and  also  the  ensuing  vortex 
sheet,  the  degree  of  which  varies  model  to  model.  The  prediction  by  Wilcox's  k-(o  model  (the 
bottom-right  figure)  seems  to  be  the  closest  to  the  measurement. 
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Figure  45.  Contour  of  axial  velocity  (U)  at  x/L  =  0.935  obtained  on  b  million  cells,  a)  measured; 

b)  SST  k-co  ;  e)  realizable  k-c\  d)  Wilcox’s  k-co 

Dynamic  Sinkage  and  Trim 

The  dynamic  sinkage  and  trim  computations  were  earned  out  using  ransFSDyMFOAM  on  a 
1.5  million-cell  hcxahedral  mesh  for  three  different  Reynolds  numbers.  Re  =  5.96  *  10f>,  1.19  x 
10  ,  1.75  x  10  '  and  the  corresponding  Froude  numbers,  Fn  =  0.138,  0.28,  0.41,  respectively. 
Unsteady  RANSE  was  solved  along  with  the  equations  of  2-DOF  (heave  and  pitch)  ship  motion. 
The  transient  runs  were  continued  until  the  resistance,  the  tnm  angle,  and  the  sinkage  reach 
statistically  steady  states.  The  time-averaged  resistance,  trim,  and  sinkage  are  shown  in  Figure 
46  along  with  the  experimental  data.  Despite  the  relatively  coarse  mesh  used  in  this  study,  the 
predictions  arc  in  fair  agreement  with  the  measurements.  More  details  on  these  calculations  can 
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be  found  in  Kim  ct  al."  . 
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Figure  46.  Prediction  of  a)  resistance,  b)  trim  and  e)  sinkage  for  the  DTMB  5415  model 


Joint  High  Speed  Sealift  (JHSS) 

This  section  involves  using  NavyFOAM’s  multiphase,  incompressible  Reynolds 
Averaged  Navier-Stokes  (RANS)  solver  for  the  Joint  High  Speed  Sealift  (JHSS)  concept  surfaee 
ship.  The  JHSS  concept  vessel  is  a  challenging  computational  ease  because  of  its  complex 
geometry  and  free  surface  interactions  w  ith  both  the  bowr  and  wraterjcts.  The  RANS  equations  are 
solved  using  NavyFOAM’s  k-omega  SST  turbulence  model.  Only  half  the  body  is  solved,  as 
symmetry  is  assumed,  and  the  domain  is  non-dimensionalized  by  ship  length.  The  coneept  vessel 
is  sealed  from  full  scale  by  constant  Froude  number  to  a  model  scale.  The  Froude  number  (Fr) 
ranges  from  -0. 23-0. 40.  The  air  and  water  phases  are  accounted  for  using  NavyFOAM’s  implicit 
Volume-of-Fluid  (VOF)  capability  The  boundary  conditions  used  for  these  computations  are: 

•  Defined  fixed  turbulent  quantities  (k,  omega,  nuTilda)  and  velocity  (U),  and  calm 
water  volume  fraction  conditions  (gamma)  at  inlet 

•  Defined  pressure  (p)  at  outlet 

•  Zero  gradient  for  all  quantities,  and  calm  water  volume  fraction  conditions 
(gamma)  at  far  field  boundaries 

•  nuTilda  =  0,  k  ~  0,  and  omega  set  to  zero  gradient  at  the  walls 

Side  and  stem  views  of  the  JHSS  geometry  can  be  seen  below  in  Figure  47.  The  left  side 
of  Figure  47  shows  the  profile  of  the  concept  vessel  including  the  gooseneck  bow,  which  has 
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very  little  clearance  from  the  free  surface.  The  right  side  of  the  figure  shows  the  waterjets  from 
the  stem  Both  the  gooseneck  bow  and  the  wateijet  inlets/nozzles  are  complex  geometry  features 
that  require  special  care/treatmcnt  in  the  meshing  process. 


Figure  47.  JHSS  concept  vessel  geometry 


Unpowered  Bare  Hull  Computations 

Initially,  computations  are  done  with  the  bare  hull  (no  wateijets  included  in  the  model)  to 
test  NavyFOAM’s  multiphase  sinkage  and  trim  capability.  Later  in  the  report  powered 
computations  involving  the  wateijets  will  be  discussed.  Figure  48  shows  the  structured  surface 
mesh  on  the  bare  hull  used  for  these  calculations.  The  surface  mesh  on  the  JHSS  displays  mesh 
refinement  around  the  free  surface  to  capture  free  surface  disturbances.  The  meshes  used  for  this 
study  were  typically  on  the  order  of  2-4M  cells  total. 


Figure  48.  JHSS  structured  surface  mesh  on  the  bare  hull 


The  computational  domain  is  split  into  many  domains  to  allow  the  computations  to  be  run 
in  parallel.  The  domain  is  split  using  the  METIS  domain  decomposition  method.  Typical  steady 
state  run  times  for  this  geometry  are  48-72  hours  depending  on  the  number  of  domain  partitions. 
Steady  state  convergence  is  assumed  when  the  forces  (pressure  and  viscous)  on  the  body  and 
sinkage  and  trim  values  change  by  a  negligible  amount  from  one  iteration  to  the  next.  Figure  49 
shows  wave  profile  results  from  NavyFOAM  computations. 
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Figure  49.  JHSS  wave  profile  on  the  hull  for  various  NavyFOAM  meshes  and  experimental 
measurements,  scaled  up  to  full  scale  (full  scale  LBP  -950  ft) 

Figure  49  shows  results  from  a  NavyFOAM  grid  study,  and  experimental  measurements 
for  a  slightly  different  JHSS  model.  A  grid  resolution  study  and  a  blocking  study  resulted  in  3 
different  meshes  (OF-A,  OF-B,  and  Ref  OF-A).  The  gnd  study  showed  that  NavyFOAM  gave 
relatively  consistent  results  for  all  three  meshes,  thus  the  grid  scheme  shown  in  Figure  48  is  used 
for  all  results  in  this  report.  Figure  49  also  shows  that  the  NavyFOAM  predicted  wave  profiles 
match  experimental  measurements  well.  However,  one  can  see  that  the  bow  wave  is  slightly 
under  predicted.  This  is  most  likely  due  to  a  lack  of  grid  resolution  in  the  bow  area  and/or  need 
for  a  sharper  volume  fraction  discretization  scheme.  Nevertheless  the  differences  are  relatively 
small  considering  the  ship  length  is  950  feet  and  the  bow  wave  is  under  predicted  by  -1  foot. 

Figure  50  shows  some  axial  velocity  boundary  layer  profiles  upstream  of  where  the 
wateijet  inlet  would  be  located  (inboard  on  the  left  and  outboard  on  the  right).  NavyFOAM 
results  (OF)  are  compared  to  both  experimental  measurements  (Exp)  and  previous  double  body 
calculations  with  another  RANS  solver  (TENASI).  The  NavyFOAM  boundary'  layers  match 
experimental  measurements  very  well.  These  boundary'  layer  plots  are  important  because 
powering  predictions  ultimately  depend  on  (amongst  other  things)  accurate  prediction  of  the 
flowficld  upstream  of  the  wateijet  inlets. 
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Figure  50.  Inboard  (left)  and  outboard  (right)  axial  velocity  boundary  layer  plots  for 
NavyFOAM  free  surface  computations  (OF),  TENASI  double-body  computations  (TEN),  and 

experimental  measurements  (Exp) 

Figure  51  shows  sinkage  and  trim  values  over  the  course  of  a  NavyFOAM  run  for  three 
Froudc  number  eases.  Each  plot  shows  a  consistent  pattern  for  both  sinkage  and  trim  for  all  run 
times,  thus  showing  that  the  multiphase  solver  is  relatively  robust,  and  wild  swings  in  sinkage 
and  trim  over  the  course  of  a  run  are  not  predicted. 
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Figure  51 .  Sinkage  and  trim  run  time  values  plotted  for  three  different  Froude  numbers 


Figure  52  shows  free  surfaee  plots  for  the  bare  hull  configuration  under  fixed  (to  the 
design  point)  and  free  sinkage  and  trim  cases.  These  plots  show  that  the  predicted  wave  profiles 
for  the  free  to  sink  and  trim  case  are  similar  to  the  profiles  predicted  under  fixed  conditions. 


Figure  52.  Fixed  and  free  sinkage  and  trim  NavyFOAM  free  surfaee  plots  colored  by  wave 

elevation 
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Figure  53  shows  total  resistance  predictions  on  the  body  over  various  Froudc  numbers  for 
both  NavyFOAM  and  experimental  measurements.  One  can  see  that  NavyFOAM  predicts  drag 
on  the  body  extremely  well  for  all  Froude  numbers  as  compared  to  experiment.  These  successful 
predictions  are  important  because  resistance  prediction  is  a  key  to  powering  predictions. 


Figure  53.  JHSS  resistance  for  various  Froude  numbers  predieted  by  experiment  and 

NavyFOAM 

Figure  54  shows  sinkage  (top)  and  trim  (bottom)  comparisons  to  experimental 
measurements  for  various  Froude  numbers.  NavyFOAM  sinkage  predictions  mateh  experimental 
measurements  very  well,  with  a  slight  diserepaney  at  the  highest  Froude  number.  NavyFOAM 
trim  predictions  also  mateh  experimental  measurements  well  over  the  range  of  Froude  numbers, 
with  a  slight  diserepaney  at  the  two  highest  Froude  numbers.  Although  the  trim  differences  may 
look  large  they  only  differ  by  fractions  of  a  degree.  In  both  the  sinkage  and  trim  eases  the  overall 
trends  are  predicted  correctly. 


Figure  54.  JHSS  bare  hull  sinkage  (top)  and  trim  (bottom)  predictions  for  various  Froude 

numbers 
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Powered  Computations  with  Waterjets 

The  rest  of  this  section  will  discuss  powered  JHSS  conditions.  The  powering  model 
described  below  is  for  the  JHSS  concept  vessel  fixed  at  design  sinkage  and  trim.  The  watcijet 
pumps  are  simulated  by  a  body  force  model  in  NavyFOAM,  which  imposes  a  pressure  jump  over 
a  volume  region  as  specified  by  the  user.  Figure  55  shows  the  surface  mesh  on  the  waterjet 
region  of  the  JHSS.  The  complexity  of  the  watcijet  inlet  geometry  led  to  the  desire  to  use 
unstructured  elements  (tetrahedral  and  prism)  to  grid  around  the  wateijet  inlets.  NavyFOAM \s 
Generalized  Grid  Interface  (GGI)  capability  allows  us  to  combine  the  structured  mesh  used  for 
the  bare  hull  computations  with  an  unstructured  grid  around  the  wateijet  inlets.  Figure  55  shows 
the  different  structured  and  unstructured  surfaces  on  the  hull. 


Figure  55.  Surface  mesh  at  the  stem  showing  GGI  region  around  wateijets 


Powered  computations  are  carried  out  similar  to  the  process  described  above  for  the  bare 
hull  case  (domain  decomposition,  steady  state  criterion,  etc.).  One  additional  steady  state 
criterion  is  added  for  the  powering  case:  the  resistance  on  the  body  must  match  an  artificial  tow 
force  (to  match  experimental  tow  tank  results)  and  the  thrust  in  the  wateijets  provided  by  the 
body  force  model  to  simulate  self  propulsion.  Figure  56  shows  axial  velocity  contours  through  an 
inboard  wateijet  inlet.  The  right  side  of  Figure  56  shows  the  volume  mesh  at  this  cross  cut.  One 
can  see  that  the  flow  remains  smooth  through  the  GGI  interfaces,  thus  NavyFOAM  s  GGI 
capability  successfully  handles  these  complex  flow  interfaces. 


Figure  56.  Axial  velocity  contours  through  the  GGI  modeled  wateijet  without  (left)  and  with 

(right)  volume  mesh  overlayed 


Figure  57  shows  axial  velocity  (Kv)  contours  predieted  by  NavyFOAM  for  the  inboard 
(right)  and  outboard  (left)  wateijets  just  upstream  of  where  the  watcijet  pumps  (or  body  force 
model  in  Navy FOAM's  case)  would  reside.  The  qualitative  look  of  the  flowfield  upstream  of  the 
pump  matches  previously  validated  TENASI  computations  very  well.  It  is  important  to 
accurately  capture  the  flow  upstream  from  the  pump  as  it  is  a  key  factor  in  final  power  (DHP) 
predictions.  NavyFOAM  predicted  self  propulsion  at  a  thrust  of  163  N,  while  experimental  tests 
resulted  in  a  self  propulsion  point  of  153  N.  There  is  a  slight  discrepancy  between  experiment 
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and  NavyFOAM  self  propulsion  points,  but  they  are  still  quite  close  considering  the  different 
methods  (tow  tank  experiments  vs.  computations)  used  in  obtaining  the  results. 


Figure  57.  Axial  velocity  contours  inside  the  wateijets 


Figure  58  shows  the  free  surface  eolored  by  wave  height  for  the  powered  JHSS  as 
predicted  by  NavyFOAM.  The  wave  pattern  is  similar  to  that  found  in  the  bare  hull  case  except 
in  the  stem  region,  as  it  should  be.  The  rooster-tail  at  the  stem  is  captured  correctly,  and  the 
affect  of  the  flow  exiting  the  wateijets  can  be  seen. 


Figure  58.  Powered  JHSS  free  surface  plot  colored  by  wave  elevation 

Figure  59  shows  a  photograph  of  the  flow  exiting  the  wateijets  during  tow  tank  tests 
compared  to  NavyFOAM’s  post-processed  results.  One  can  see  that  the  NavyFOAM 
computations  predict  the  complex  physics  taking  place  at  the  stem  very  well.  The  interaction 
between  the  rooster-tail  shooting  out  from  underneath  the  stem  and  the  flow  exiting  the  wateijets 
is  captured  very  well  in  NavyFOAM. 
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Figure  59.  Experimental  photograph  (left)  and  NavyFOAM  post-processed  JHSS  powered  stem 

In  conclusion,  NavyFOAM  displays  the  capability  to  accurately  predict  complex  physics 
for  surface  ships.  The  JHSS  concept  vessel  is  an  especially  complex  surface  ship  due  to  the 
gooseneck  bow  and  the  wateijet  inlets/powering.  Initially,  the  hull  resistance,  sinkage  and  trim  of 
the  bare  hull  case  are  predicted  with  the  hull  free  to  sink  and  trim.  Navy  FOAM  results  match 
experimental  measurements  very  well  for  various  Froude  numbers  that  take  the  concept  vessel 
through  different  ship  attitudes.  Finally,  self  propulsion  is  predicted  with  a  body  force  model  (in 
place  of  the  actual  waterjet  pumps)  providing  thrust  that  balances  out  the  model's  resistance. 
NavyFOAM  powering  predictions  match  experimental  measurements  well. 

Summary 

This  report  provides  a  guide  to  NavyFOAM  VI. 0,  which  is  based  on  the  OpcnFOAM 
open  source  software.  A  brief  technical  description  of  the  code  is  given  w  ith  an  emphasis  on 
those  changes  made  for  NavyFOAM  V1.0  that  differentiates  it  from  the  standard  OpenFOAM 
offering.  More  details  on  OpenFOAM  specifically  can  be  found  in  the  OpenFOAM  guides 
referenced  in  this  report.  Complementing  the  technical  description  of  NavyFOAM  changes 
there  is  a  User’s  Guide  section  to  help  users  of  NavyFOAM  properly  implement  the  use  of  these 
improvements.  Results  are  presented  for  a  variety  of  Navy  relevant  configurations  including  a 
fully  submerged  axisymmetric  body,  a  tanker,  DTMB  Model  5415  (pre-contract  design  for  the 
DDG-51)  and  the  Joint  High  Speed  Sealift  (JHSS)  concept  Results  for  all  cases  compare  will 
with  experimental  data.  Results  are  also  presented  for  a  variety  of  grid  types  and  turbulence 
models  providing  some  indication  of  the  capabilities  available  with  the  code.  Finally,  several 
tutorials  as  well  as  other  information  that  is  directly  aimed  at  helping  users  successfully  use  the 
code  are  also  provided. 
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Appendix  A:  Supplemental  User’s  Guide 

This  section  is  meant  to  supplement  the  original  OpenFOAM  User’s  Guide.  For  more 
detailed  information  on  the  OpenFOAM  code  and  settings  consult  the  User’s  Guide: 
http://foam.soureeforge.net/doe/Guides-a4/UserGuide.pdf. 

Standard  solvers  that  are  used: 

•  ieoFoam  (transient,  laminar,  incompressible,  single  phase) 

•  simpleFoam  (steady,  RANS,  incompressible,  single  phase) 

•  interFoam  (transient,  RANS,  incompressible,  multi-phase) 

Files  Contained  in  the  system  Directory  include: 

(1)  eontrolDict 

(2)  deeomposeParDiet 

(3)  fvSehemes 

(4)  fvSolution 

eontrolDict 

Dictionary  that  controls  run  parameters  and  output  data  from  the  run. 

Application:  here  you  specify  the  solver  used  (i.e.  simpleFoam,  ieoFoam,  etc.).  This  is 
not  essential.  Specifying  the  executable  here  doesn’t  seem  to  have  any  control  over  the  run. 

Start  From:  here  you  specify  when  to  start  the  run  from.  Options  are: 

firstTime  -  start  from  the  earliest  time  directory  available 

startTime  -  start  from  the  time  specified  by  startTime  on  the  next  line 

latestTime  start  from  the  most  reeent  time  directory  available. 

startTime :  enter  the  time  to  start  the  run  from  when  startFrom  ?  startTime;  is  selected 

StopAt :  here  you  specify  when  to  stop  the  application.  Options  are: 

endTime  -  stop  run  at  time  specified  by  endTime  on  the  next  line 

writeNow  -  stop  the  run  at  the  next  iteration  and  write  out  the  last  time  step 

noWriteNow  -  stop  the  run  at  the  next  iteration  and  don't  write  time  info  at  last  time 

step 

nextWnte  -  stop  the  run  at  the  next  scheduled  write  time  specified  by  writeControl 

endTime:  time  when  run  will  stop.  Only  valid  when  stopAt  ?  endtime;  is  selected 

deltaT:  here  you  specify  the  run  time  step.  This  is  typieally  /  for  steady  runs  and  for 
unsteady  runs  this  is  superceded  by  maxCo ,  which  is  discussed  below. 

WriteControl:  here  you  specify  when/how  to  output  information.  Options  are: 

timeStep  writes  data  every  writelnterval  time  steps 

runTime  writes  data  every  writelnterval  seconds  of  run  time 
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adjustableRunTime  -  this  is  the  same  as  runTime,  exeept  it  will  adjust  the  time  steps  to 
coincide  with  the  writelnterval 

cpuTime  -  writes  data  every  writelnterval  seconds  of  CPU  time 

clockTime  -  w  rites  data  every  writelnterval  seconds  of  real  time 

writelnterval :  scalar  data  that  specifies  time  output  from  writeControl 

purgeWrite :  this  is  an  integer  value  that  specifies  how  many  output  time  directories  arc 
stored.  When  the  max  number  of  purgeWrite  directories  have  been  written,  the  newest  time  data 
will  over  write  the  oldest  time  directory.  A  value  of  0  means  that  no  time  data  will  be 
overwritten. 

WriteFormat :  here  you  specify  the  format  of  the  output  data.  Options  arc: 
ascii  -  ASCII  data  with  the  amount  of  significant  figures  specified  by  write  Precision 
binary  -  Binary  data 
ascii  is  typically  used 

writePrecision :  number  of  significant  figures  ASCII  w  rite  Format  is  written  out  to. 
WriteCompression :  here  you  specify  the  compression  (if  any)  of  the  output  files.  Options 

are: 

uncompressed  No  data  compression 
compressed  -  gzip  compression 
eompressed  has  typically  been  used 

timeFormat:  here  you  speeify  the  format  for  the  time  directory  names.  Options  are: 
fixed  -  all  time  directories  are  written  in  fixed  format  (i.e.  123.456) 
seientifie  -  all  time  directories  are  labeled  in  scientific  format  (i.e.  1 .23456e+03) 
general  -  specifies  seientifie  format  of  exponent  is  less  than  -4,  otherwise  fixed 
general  has  typically  been  used 

timePrecision :  here  you  speeify  the  number  of  significant  figures  for  timeFormat  time 
directory  labels.  6  is  the  default  and  is  typically  used 

GraphFormat :  here  you  speeify  the  graph  data  written  by  an  application.  Options  arc: 

raw  -  data  in  ASCII  format  in  columns 

gnuplot  -  data  in  gnuplot  format 

xmgr  -  data  in  Grace/xmgr  format 

jplot  -  data  in  jPlot  format 

typieall  raw  is  specified,  but  this  depends  on  the  user  preference. 

RunTime  Modifiable :  here  you  specify  yes  to  have  the  directories  (ineluding  controlDict) 
read  at  the  beginning  of  each  time  step,  or  no  to  have  the  run  proceed  without  rereading 
directories.  This  is  typically  set  to  yes,  and  should  remain  so  to  avoid  eonfusion. 
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AdjustTimeStep :  select  yes  for  an  adjustable  time  step  (whose  size  is  determined  by 
maxCo  and  maxDeltaT ),  or  no  for  a  constant  time  step  as  specified  by  deltaT. 

MaxCo:  This  adjusts  the  time  step  to  aehieve  the  specified  maximum  CFL  number.  This 
value  is  typically  low  0(1)  at  the  beginning  of  a  run,  and  ramped  up  0(10)  once  the  solution 
becomes  more  stable. 

MaxDeltaT :  this  value  is  a  limit  to  the  maximum  time  step  achievable,  when  a  maxCo  is 
specified. 

**  At  the  bottom  of  the  controlDict  (or  anywhere  within)  you  can  speeify  additional 
libraries  or  functions  to  be  loaded  at  run-time.  Turbulence  model  and  dynamic  mesh  motion  libs , 
and  force  functions  are  common  examples  of  what  have  been  used. 

decomposeParDict 

Dictionary  that  contains  all  input  information  on  domain  decomposition  for  parallel 
processing  runs. 

Subdietionary:  numberOfDomains 

Here  you  specify  the  number  of  partitions  you  want  your  grid  to  get  split  up  into. 

Subdietionary:  method 

Here  you  specify  the  method  of  domain  decomposition.  Options  are: 

simple  domain  is  decomposed  into  volumes  that  are  similar  in  all  coordinate 
directions. 

hierarchical  -  user  can  speeify  the  order  of  direction  for  simple  decomposition. 

scotch  -  attempts  to  minimize  the  number  of  geometric  boundaries. 

metis  similar  to  scotch,  but  is  not  free  for  commercial  use,  it  will  eventually  be 
discontinued. 

manual  -  user  manually  specifies  decomposition  of  eaeh  cell  to  a  processor. 

metis  has  traditionally  been  used  during  testing,  occasionally  simple  decomposition  is 
used.  User  should  use  metis  except  for  special  eireumstanees. 

Subdietionary:  simpleCoeffs 

Here  you  speeify  n  number  of  processors  in  eaeh  direction  to  decompose  domain,  and 
cell  skew  factor  (delta)  for  decomposition.  This  is  only  necessary  for  simple  domain 
decomposition. 

Subdietionary:  hierarehiealCoeffs 

Here  you  speeify  n  number  of  processors  in  eaeh  direction  to  decompose  domain,  and 
cell  skew  factor  (delta),  and  the  order  of  directions  (i.e.  xyz  or  zyx)  for  decomposition.  This  is 
only  necessary  for  hierarchical  domain  decomposition. 

Subditionary:  scotchCoeffs 

Here  you  speeify  the  weighting  factors  ( processorlVeights )  for  eaeh  individual  processor. 
The  numbers  for  eaeh  processor  are  normalized,  so  any  values  can  be  accepted,  no  matter  the 
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sum.  You  can  also  specify  the  strategy,  but  it  is  not  clear  what  this  value  means.  This  is  only 
necessary  for  scotch  domain  decomposition. 

Subdictionary:  metisCoeffs 

Here  you  specify  processorWeights ,  as  described  above  in  scotchCoeffs.  This  is  only 
necessary  if  you  specify  metis  decomposition. 

Subdictionary:  manualCoeffs 

Here  you  specify  the  name  of  a  data  file  that  contains  processor  allocations  for  each  cell. 
This  is  only  necessary  if  you  specify  manual  decomposition. 

Subdictionary:  distributed 

This  is  an  optional  subdictionary,  where  you  state  yes  or  no,  whether  there  is  geometry 
data  in  any  other  directories  that  needs  to  be  decomposed  with  the  current  directory. 

Subdictionary:  roots 

If  you  chose  yes  to  distributed ,  here  you  list  the  address(cs)  to  additional  directories. 

fvSchemes 

Dictionary  that  contains  numerical  scheme  input:  interpolation  methods,  temporal  and 
spatial  discretization  information 

Subditctionary:  ddtSchemes 

This  subdictionary  contains  first  time  derivative  discretization  method 

Eider  (lsl  order  implicit)  is  the  only  method  that  has  been  consistently  used 

Subdictionary:  gradSchemes 

This  subdictionary  contains  discretization  information  for  the  gradient  terms 

Gauss  linear  and  leastSquares .  The  Gauss  linear  method  has  been  used  most  frequently, 
yielding  consistent  results.  The  leastSquares  method  is  believed  to  be  more  accurate  when 
calculating  the  gradient  on  non-uniform  meshes,  but  bugs  were  encountered  early  in  the  V  &  V 
process.  Improvement  of  the  leastSquares  gradient  method  has  been  an  ongoing  effort. 

Limiting  is  available  (i.e.  Gauss  linear  limited );  however,  this  was  proven  to  negatively 
affect  the  solution  during  the  V  &  V  process. 

Subdictionary;  divSchemes 

This  subdictionary  contains  discretization  information  on  the  divergence  terms. 

Div(phi,U)  is  commonly  refered  to  as  the  convection  term  in  the  momentum  equation 
Typically,  a  2nd  order  upwind  is  used  for  this  term,(7a//^  linearUpwind  cellLimited  Gauss  linear 
1.0. 

Turbulent  quantities  (nuTilda,  k,  omega,  epsilon,  etc.)  are  usually  solved  1st  order  upwind 
(Gauss  upwind ),  but  it  may  be  necessary  to  solve  the  2nd  order  upwind  (Gauss  linearUpwind 
cellLimited  Gauss  linear  1.0). 

Multiphase  gamma  terms. . . 
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*Note:  There  will  be  an  error  message  if  the  term  div((nuEjf*dev(grad(U).  T())))  is  left 
out,  or  if  a  divScheme  is  applied  to  this  term.  Instead  this  term  needs  to  be  included  with  a 
gradient  scheme  discretization.  Gauss  linear  is  the  most  common  scheme  used  with  this  term. 

Subdictionary:  laplacianSchemes 

This  subdictionary  contains  discretization  information  for  the  laplaeian  terms. 

Gauss  linear  corrected  has  been  the  standard  scheme  used.  Limiting  ean  be  used  ( Gauss 
linear  limited  0.0  =  uncorrectcd  and  Gauss  linear  limited  1.0  =  eorreeted).  The  correction  refers 
to  treatment  of  non-orthogonal  terms.  An  uneorreeted  solution  (Gauss  linear  limited  X.X  with 
X.X <  1  will  not  converge  to  the  same,  more  aeeurate,  solution  as  Gauss  linear  corrected). 

Subdictionary:  interpolationSchemes 

This  subdietionary  contains  information  on  interpolation  that  is  usually  from  cell  center 
to  cell  face. 

Linear  has  been  the  standard  scheme  used;  however,  reconCentral  seheme  is  being 
developed  and  is  believed  to  be  more  accurate  for  meshes  with  significant  non-orthogonality 
(unstructured  meshes). 

Subdietionary:  snGradSehemes 

This  subdietionary  contains  information  on  surfaee  normal  gradient  terms.  This  term 
specifies  the  portion  of  the  gradient  at  a  cell’s  face  that  is  normal  to  the  faee. 

corrected  has  been  the  standard  seheme  used  for  this  subdictionary  and  other  schemes 
have  not  been  used  significantly. 

Subdietionary:  fluxRequired 

This  subdictionary  contains  information  for  variables  whose  flux  is  calculated  in  the 
application. 

The  llux  is  required  for  pressure  (p  and  pd)  for  most  calculations,  beeause  a  pressure 
equation  is  solved.  Thus  the  default  is  usually  set  to  no  and  the  variable  p  or  pd  is  specified  in 
the  subdietionary,  so  that  the  flux  is  ealeulated  for  the  pressure. 

fvSolution 

Dietionary  that  contains  algorithm  and  linear  system  solvers  information,  sueh  as  solver 
settings  and  toleranees  for  convergence.  The  segregated  solver  variables  and  solution  algorithms 
will  vary  with  ehoiee  of  problem  solver  executable  (i.e.  simpleFoam  and  interFoam  require 
different  settings). 

solvers  contains  the  linear  system  of  equation  algorithms  and  tolerances  for  all  the 
variables  (OpenFoam  uses  a  segregated  solver). 

The  common  linear  solvers  tested  and  used  are  conjugate  gradient  solvers  (PCG/PBiCG) 
and  multi-grid  solvers  (GAMG/AAMG). 

Many  tolerance  settings  have  been  used,  but  eommon  values  tested  have  been: 

P  -  tolerance  1  e- 1 0;  relTol  0.0 1 ;  and  minlter  =  1 ; 

U,  gamma,  and  turbulent  quantities  -  toleranee  le-07;  rcITol  0.0;  and  minlter  =  1 ; 
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The  tolerance  criterion  is  satisfied  when  the  residual  for  the  linear  solver  reaches  the 
prescribed  value.  The  rclTol  criterion  is  satisfied  when  the  linear  solver  residual  has  dropped  by 
the  order  of  magnitude  specified  by:  1.0/relTol  (i.e.  relTol  =  0.01  corresponds  to  the  residual 
dropping  two  orders  of  magnitude  or  10~). 

Note:  be  sure  to  use  the  mini  ter— I;  command  for  all  the  linear  solvers,  as  some  solution 
tolerances  may  be  set  such  that  variables  will  stop  iterating  prematurely  in  the  solution  process, 
thus  leading  to  inaccurate  solutions  that  may  look  OK, 

Pressure  algorithms  (SIMPLE  and  PISO)  - 

For  simpleFoam  (steady  state,  single  phase,  RANS  solver)  the  solver  variables  needed  are 
p,  U,  turbulent  quantities  (nuTilda  for  SA,  k  and  omega  for  SST  or  WKO,  k  and  epsilon  for  k- 
epsilon,  etc.). 

The  pressure  algorithm  should  be  SIMPLE,  which  includes: 

nNonOrthogonalCorrectors  #;  Where  the  #  selected  will  determine  the  number  of 
additional  pressure  solver  loops.  For  example,  for  a  value  of  #  =  1,  the  solver  will  iterate  over  the 
pressure  equation  twice.  For  runs  with  meshes  of  good  quality  these  additional  loops  are  not 
needed.  However,  for  meshes  containing  a  lot  of  skew  or  nonOrthogonality,  values  between  1 
and  3  will  add  stability  (as  well  as  additional  computational  time)  to  the  solution. 

nCorrectors? 

pRefCell  and  pReJValue  must  both  be  set  or  errors  will  occur.  These  values  represent 
what  the  reference  pressure  value  (most  likely  frec-strcam)  and  a  cell  number  where  this  value 
occurs  (a  cell  located  in  the  free  stream).  All  previous  runs  have  used  values  of  0  for  both 
without  problem. 

For  unsteady  and  pseudo  time-marching  steady  solvers  like  interFoam  and 
ransNavylnterFoam 

The  pressure  algorithm  should  be  PISO,  which  includes: 
nNonOrthogonalCorreetors 

relaxationFactors  contains  the  under-relaxation  values  for  all  the  linear  solver  variables 
(U,p,pd,omega,k, gamma,  etc.).  The  relaxation  factors  correspond  to  : 

0.0=  Fully  relaxed 

0.0  <  #  <  1.0  =  variable  under-relaxation 

1 .0  =  No  relaxation 

For  all  implicit  and  most  explicit  runs,  under-relaxation  should  be  used  on  all  variables  to 
varying  degrees.  The  pressure  terms  (p  for  simpleFoam,  pd  for  other  solvers)  require  more 
undcr-relaxation  than  other  variables  and  this  value  is  usually  very  important  to  solution  stability. 
It  is  not  uncommon  to  start  a  calculation  out  with  a  value  of  0.1  or  0.2  and  wait  for  the  initial 
solution  to  develop  and  then  ramp  the  value  up  progressively  to  0.3.  The  velocity  (U)  usually 
starts  around  0.5  and  ramps  up  to  0.7  or  0.8  as  the  solution  stabilizes.  The  turbulent  quantities 
and  gamma  usually  affect  the  solution  startup  less  and  typical  values  are  0.6  and  0.7-0. 8 
respectively. 
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Appendix  B:  Utility  Programs 

Several  programs  have  been  written  to  simplify  post-processing  NavyFOAM  output  data. 
Descriptions  of  their  use  arc  given  here. 

dataFunkyFieldComparison 

Description 

This  utility  reads  the  scalar  and  vector  volume  field  from  FOAM  numerical  solution 
data  files  and  compares  them  with  the  analytic  solutions  specified  by  a  dictionary  file  named 
funkyFieldsDict  stored  under  directory  $ CA SE_ DIR/system.  The  utility  will  first  match  the  field 
names  specified  in  funkyFieldsDict  with  those  stored  in  the  solution  files  and  then  calculate  the/^ 
and  L2  norm  of  the  error  and  send  the  report  to  standard  output  (stdout).  The  error  will  only  be 
calculated  if  the  field  name  is  found  in  both  funkyFieldsDict  and  solution  files. 

This  utility  is  part  of  NavyFOAM. 

Usage 

The  command  line  usage  looks  like 

dataFunkyFieldComparison  [-case  dir]  [-time  time]  [-latestTime] 

The  optional  options  arc 
-case  dir:  specifics  the  case  directory; 

-time  time:  selects  the  time  step; 

-latestTime:  selects  the  latest  time  step. 

Without  any  option,  the  utility  will  read  the  data  files  for  all  time  steps  stored  under  the 
current  case  directory  . 

Installation 

1)  Create  a  working  copy  using  svn  checkout.  The  recommended  local  directory  to 
checkout  the  package  is 

NavyFOAM/apphcations/utihties/postProcessing 

svn  checkout  svn+ssh://blackwater/5700/NavyFOAM/IniegrationBranches/NavyF()AM-1.5-dcv- 
rev995/NavyFOAM/applications/utilities/postProcessing/dataComparison 

2)  Go  to  directory  JdataComparison/dataFunkyFieldComparison  and  compile  the  package: 

wmake 

The  compiler  may  generate  some  warning  message  which  can  be  ignored  in  this  case. 
The  generated  executable  file  dataFunkyFieldComparison  can  be  found  in  a  user  application  binary 
file  directory  specified  by  $FOAM_ USE R_APPB/N. 

Warning:  To  compile  this  utility  at  least  version  2.1  of  Bison  has  to  be  installed.  Check  w  ith 
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bison  -V 


on  the  command  line  before  trying  to  compile  it.  Go  to  httn://www.unu.oru/softwarc/bisoa  for 
more  information  regarding  Bison. 

Expression  Syntax 

The  example  below  shows  some  useful  expression  syntax.  The  most  complete 
documentation  of  the  expression  syntax  is  the  source  file  for  the  Bison-grammar  in  the  package: 
ValueExpressionParser.yy 
ValueExpressionLexer.  // 

Example 

The  dictionary  file  funky FieldsDict  used  in  a  Taylor  vortex  test  case  is  shown  below: 

FoamFile 

{ 

version  2.0: 
format  ascii; 
class  dictionary; 
object  funkyFieldsDict; 

* 

t 

// *************************************// 

expressions 

( 

TaylorVortex  Velocity 

i 

field  U; 

expression  "vector(sin(pos().x)*cos(pos().y),  -cos(pos().x)*.sin(pos().y),  0)*exp(-0.2*time())"; 

i, 

TaylorVortexPressure 

{ 

field  p; 

expression  "0.25*(cos(2.*pos().x)+cos(2.*pos().y))*pow(exp(-0.2*time()),2)"; 

i 

r 

); 

The  above  dictionary  file  specifies  the  analytic  solution  of  the  Taylor-Green  vortex 
problem: 

p{.\ ,  y,/)  =  —  (cos2.v  +  cos2y)F2(/) 

4 

u(x,y,l)  =  sin.v  cos yF(t) 
v(x ,yj)  =  -cosx  sinyF(/) 

with  F(t)  =  e~2"  and  v  =  0. 1 . 

Suppose  we  already  have  the  solution  files  at  time  =  1  as  the  latest  time  step  stored 
under  SCASE  DIR/1.  The  command  line 

dataFunkyFieldComparison  -latestTime 
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will  yield  the  output  shown  below: 


ummmmmmmmmmmmnmmmmmmummmmmmn 


Time  =  1 

mmMMMmmmffltMmmMMumtMMMmtmMmitM 


volScalarField;  p 


max  error  occurs  at  cell:  0 
maxErr  =  0.0004091 10343 
rmsErr :  0.000137916310931 


volVectorField:  U 


max  error  in  magnitude  occurs  at  cell:  12216 
maxErr  =  6.79861 100733e-05 
rmsErr  =  2  4499607940  le-05 

max  error  in  X-component  occurs  at  cell:  7680 
maxErr  =  3.4175343c-05 
rmsErr  =  1 .447386995 1 3e-05 


max  error  in  Y-component  occurs  at  cell:  12217 
maxErr  =  6.701 7086c-05 
rmsErr  2.1343033345e-05 


max  error  in  Z-component  occurs  at  cell:  591 1 
maxErr  =  0 
rmsErr  =  0 


For  a  scalar  volume  field  such  as  the  pressure  field,  the  maximum  error  (maxErr)  and 
root-mean-squares  error  (rmsErr)  are  calculated  as  follows: 


maxErr  =  L,  (pn  -  pn0 )  =  max  (|  p"  -  pn0.  |) 

\<,i  <*Np 


where  the  subscript  “0”  denotes  the  analytic  solution,  the  superscript  represents  the  n- 
th  time  step,  and  NP  is  the  total  number  of  cells  for  cell-centered  finite  volume  method 

For  a  vector  volume  field  such  as  the  velocity  field,  the  maxErr  and  rmsErr  are  calculated 
for  the  vector  magnitude: 


rmsErr  =  A, (|  V  |"  -|K0|")  = 


maxErr  =  Z.„(|  V  |"  -\VQ  |")  = 
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and  for  each  component: 
x-component 


maxErr  =  L„(u"  -i/„)  -  max  (|w"  -ua0j  |) 

15  I  5<V,, 


rmsErr  =  L2  («"  -  Hg )  = 


^’-component 


maxErr  =  L„(v"  -  )  =  max  (|v"  -  |) 

15  /  £Nr 


rmsErr  =  L2  (v"  -  v0" )  =  v,"  -  ) 3 


z-component 


maxErr  =  L„(wn  -w”)=  max  (K  “  woj  I) 

»  <iVr 


rmsErr  =  L2(wn  -  vv")  = 


An  optional  option  patches  will  be  added  to  calculate  the  and  L2  nonrj  of  the  error  in 
patch-intemal-field  (cells  that  directly  connecting  the  patch)  for  specified  patches. 

NavyFOAMToTecpIot 

Description 

This  utility  reads  the  scalar  and  vector  volume  and  boundary  patch  fields  from  FOAM 
numerical  solution  data  files  and  converts  them  to  Tccplot  data  files. 

This  utility  is  part  of  Navy  FOAM. 

Usage 

The  command  line  usage  looks  like 

NavyFoamToTecpIot  [-region  name]  [-case  dir]  [-fields  fieldsList] 

[-patches  patchesList]  [-time  time]  [-latestTime] 

The  optional  options  are 

-region  name:  specifies  the  region  name; 

-case  dir:  specifies  the  case  directory; 

-fields  fieldsList:  specifies  a  list  of  fields  to  output,  e.g.  ‘(  p  IJ  gamma  )\  or  4(U)\ 
-patches  patchesList:  specifics  a  list  of  patches  to  output,  c.g.  ‘(inlet  outlet  wall)’,  or 


‘(hull)’; 


-time  time:  selects  the  time  step; 

-latestTime:  selects  the  latest  time  step. 

Without  any  option,  the  utility  will  read  the  data  files  for  all  time  steps  stored  under 
the  current  case  directory  and  convert  the  volume  field  to  Tecplot  data. 
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Installation 


1)  Create  a  working  copy  using  svn  checkout.  The  recommended  local  directory  to 
checkout  the  package  is 

NavyFOAM/applications/utilities/postProcessing/dataConversion 

svn  checkout  svn+ssh://blackwater/5700/NavyFOAM/IntegrationBranches/NavyFOAM-l  .5-dev- 
revl745/NavyFOAM/applications/utilities/postProcessing/dataCon  version 
/NavyFoamToTecplot 

2)  Go  to  directory  VdataConversion/NavyFoamToTecpIot  and  compile  the  package: 

wmake 

The  generated  executable  file  NavyFoamToTecplot  can  be  found  in  a  user  application 
binary  file  directory  specified  by  $FOAMJJSER_APPBIN. 

Output 

The  output  files  can  be  found  in  $CASE_DIR/TecplotData,  for  example 
NavyFoamTo'I ecplot  -fields  4(p  U)’  -patches  ‘(hull)’  -time  10 

will  generate  10.dat  and  hull_10.dat  in  $CASE_DIR/TecplotData. 


NavyCellSetToTecpIot 

Description 

This  utility  reads  a  user  specified  cellSet  file  in  ./constant/polymcsh/scts/  and  convert 
it  to  Tecplot  data  file. 

This  utility  is  part  of  NavyFOAM. 

Usage 


The  command  line  usage  looks  like 

NavyCellSetToTecpIot  <cellSetFileName>  [-region  name]  [-case  dir] 
[-patches  patchesList] 


The  mandatory  argument  is 

cellSetFileNamc:  specifies  the  file  name  for  the  cellSet,  c.g.  wavcDampingCclls; 

The  optional  options  are 

-region  name:  specifics  the  region  name; 

-case  dir:  specifies  the  case  directory; 
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-patches  patchesList:  specifics  a  list  of  patches  to  output,  e.g.  ‘(  inlet  outlet  wall  )\  or 

‘(hull)’; 

Installation 

1)  Create  a  working  copy  using  svn  checkout.  The  recommended  local  directory  to 
checkout  the  package  is 

NavyFOAM/applications/utilities/postProcessing/dataConversion 

svn  checkout  svn+ssh://blackwater/5700/NavyFOAM/IntegrationBranches/NavyFOAM- 1 .5-dev- 
revl  745/Navy  FOAM/applications/utilities/postProcessing/dataConversion 
/NavyCellSetToTecpIot 

2)  Go  to  directory  ./dataConversion/NavyCellSetToTecpIot  and  compile  the  package: 

wmake 

The  generated  executable  file  NavyCellSetToTecpIot  can  be  found  in  a  user 
application  binary  file  directory  specified  by  $FOAMJJSER_APPBIN. 


Output 


The  output  files  can  be  found  in  $CASE_DIR/TecplotData,  for  example 
NavyCellSetToTecpIot  waveDampingCells  -patches  ‘(hull  symmetry)' 

may  generate  waveDampingCells.dat  in  $CASE  DIR/TecpIotData.  The  Tecplot  data  file 
waveDampingCells.dat  should  contain  the  following  zones: 

ZONE  T  =  vo  I  Mesh 
ZONE  T  =  hull 
ZONE  T  =  symmetry 
ZONE  T  waveDampingCells 


The  first  zone  is  the  volume  mesh.  The  next  two  zones  are  the  surface  mesh  for  the 
user  specified  patches.  The  last  zone  is  a  volume  mesh  for  the  cellSet.  In  Tecplot,  the  mesh 
of  each  zone  can  be  plotted  in  different  colors. 

NavyFaceSetToTecpIot 

Description 

This  utility  reads  a  user  specified  faceSet  file  in  ./constant/polymesh/scts/  and  converts  it 
to  Tecplot  data  file. 

This  utility  is  part  of  NavyFOAM. 

Usage 


The  command  line  usage  looks  like 
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NavyFaceSetToTecplot  <faceSetFileName>  [-region  name]  [-case  dir] 
[-patches  patchesList] 


The  mandatory  argument  is 

faccSetFilcName:  specifies  the  file  name  for  the  faceSet,  c.g.  skcwFaces; 

The  optional  options  arc 

-region  name;  specifies  the  region  name; 

-case  dir:  specifies  the  case  directory; 

-patches  patchesList:  specifies  a  list  of  patches  to  output,  e.g.  fc(  inlet  outlet  wall  )\  or 

‘(hull)’; 

Installation 

1)  Create  a  working  copy  using  svn  checkout.  The  recommended  local  directory  to 
checkout  the  package  is 

NavyFOAM/applications/utilities/postProcessing/dataConversion 

svn  checkout  svn+ssh://blackwater/5700/NavyFOAM/lntegrationBranches/NavyFOAM-l  .5-de\- 
revT745/NavyFOAM/applications/utilities/postProcessing/dataConversion 
/N  avy  F  aceSetToT  ecplot 

2)  Go  to  directory  VdataConversion/NavyFaceSetToTecpIot  and  compile  the  package: 

wmake 

The  generated  executable  file  NavyFaceSetToTecplot  can  be  found  in  a  user 
application  binary  file  directory  specified  by  $FOAM_USER_APPB!N. 

Output 


The  output  files  can  be  found  in  $CASE_DIR/TecpiotData,  for  example 

NavyFaceSetToTecplot  skewFaces  -patches  ‘(hull  symmetry)’ 

may  generate  skewFaces.dat  in  $CASE_DIR/TecpiotData.  The  Tecplot  data  file 
skewFaces.dat  should  contain  the  following  zones: 

ZONE  T  -  volMesh 
ZONE  T  =  hull 
ZONE  T  =  symmetry 
ZONE  T  =  skewFaces 

The  first  zone  is  the  volume  mesh.  The  next  two  zones  are  the  surface  mesh  for  the 
user  specified  patches.  The  last  zone  is  a  surface  mesh  for  the  faceSet.  In  Tecplot,  the  mesh 
of  each  zone  can  be  plotted  in  different  colors. 
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Appendix  C:  icoFoam  Lid  Driven  Cavity  Tutorial 


This  tutonal  involves  using  the  laminar,  transient,  incompressible  solver  for  a  2-D 
cavity.  The  cavity  consists  of  4  walls,  where  the  top  wall  is  moving,  and  the  other  walls  arc 
stationary  .  First,  we  will  go  over  pre-processing  and  case  setup,  then  we  will  run  the  test  case, 
and  Finally  we  will  look  at  some  post-processed  results. 

For  more  detailed  information  on  the  OpenFOAM  code  and  settings  consult  the  User's 
Guide:  http://foam.sourccforge.nct/doc/Guides-a4/UscrGuide.pdf. 

Pre-Processing  and  Case  Setup 

Upon  looking  in  the  case  directory  (screen  capture  below)  you  will  notice  that  there  are 
directories  labeled  system  and  0 ,  and  a  File  named  transportProperties.  There  is  also  a  File 
labeled  2D  cavity  allCoarseStr.cas,  which  is  a  mesh  exported  from  Gridgcn  in  Fluent  ASCII 
double  precision  format.  We  will  discuss  the  existing  Files  and  directories  later,  but  now  we  need 
to  import  the  fluent  .cas  File  into  OpenFoam.  This  will  be  done  using  the  command 
fluentMesh  To  Foam . 

[kdelaney^Retech  allCoarseStruct ] £  1 
total  1.2H 

-rw-rw-r--  1  kdelaney  kdelaney  1.2M  Feb  19  12:31  2D_cavity_allCoarseStr .cas 

drwxrwxr-x  2  kdelaney  kdelaney  4. OK  Mar  9  16:53  systeu 

-rw-r -  1  kdelaney  kdelaney  886  Apr  1  08:35  transportProperties 

drwxr-xr-x  2  kdelaney  kdelaney  4. OK  Apr  1  08:35  0 

Mesh  Input 

Now  enter  the  JluentMeshToFoam  command  as  seen  below  and  view  the  output  from 
the  screen  dump.  The  extra  “|  tee  conversionLog"  is  not  essential  and  is  only  included  to  record 
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the  screen  dump  in  a  file  named  conversionLog.  Your  output  should  be  the  same  as  the  screen 
captures. 

There  is  a  lot  of  information  on  the  screen  dump,  most  of  which  is  self  explanatory.  The 
most  important  part  to  notice  is  the  last  two  lines  of  text  in  the  screen  dump  which  tell  you  that 
the  mesh  information  has  been  written  into  a  directory  named  polyMesIt  inside  a  newly  created 
directory  named  constant.  Finally  the  command  ends  successfully  with  the  “End.”  statement. 
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[kdelaneyiSRetech  allCoarseStruct 1$  f luentMeshToFoar  2D_cavity_allCoarseStr . cas  |  tee  conversionLog 

/* - *\ 


w 

/ 

F 

ield 

|  OpenFOAM: 

The  Open  Source  CFD  Toolbox 

w 

/ 

0 

peration 

|  Version: 

1 . 5-dev 

\\ 

/ 

A 

nd 

|  Web: 

http: //www . OpenFOAM . org 

w/ 

M 

anipulat ion 

Exec  :  f luentMeshToFoan  2D_cavity_allCoarseStr . cas 

Date  :  Mar  09  2010 

Tine  :  16:40:45 

Host  :  Retech 

P1D  :  5119 

Case  :  /hotie/kdelaney/NavyFoan_runs/ run/f lat_plate/rovingWal ICavity/a 1 ICoarseStruc t 

nProcs  :  1 

//..**...*♦**♦.♦**********.*  .  ******  // 

Create  tire 

Reading  header:  “exported  fror  Cridgen  15.14R1" 

Dirension  of  grid:  3 
Nurber  of  points:  9522 
nurber  of  faces:  18632 
Nurber  of  cells:  4624 
Reading  points 

Other  readCe HCroupData :  2  1  1210  1  4 
Reading  uniforr  cells 

Read  zonel:2  nare: fluid  patchTypelD : fluid 
Reading  zone  data 

Erbedded  blocks  in  corrent  or  unknown:  ( 

Found  end  of  section  in  unknown:) 

.Reading  uniforr  faces 

Read  zonel:3  nare : interior-3  patchTypelD : inter lor 
Reading  zone  data 

Erbedded  blocks  in  corrent  or  unknown:  ( 

Found  end  of  section  in  unknown:) 

.Reading  uniforr  faces 

Read  zonel:4  nare:lnlet-4  patchTypelD :wa 1 1 
Reading  zone  data 

Erbedded  blocks  in  corrent  or  unknown:  ( 

Found  end  of  section  in  unknown:) 

.Reading  uniforr  faces 

Read  zonel:5  nare:Top-5  patchTypelD : wall 
Reading  zone  data 

Erbedded  blocks  in  corrent  or  unknown:  ( 

Found  end  of  section  in  unknown:) 

.Reading  uniforr  faces 

Read  zonel:6  nare: side_l-6  patchTypelD :wal 1 
Reading  zone  data 

Erbedded  blocks  in  corrent  or  unknown:  ( 

Found  end  of  section  in  unknown:) 

.Reading  uniforr  faces 

Read  zonel:7  nare : Out  let-7  pa tchTypelD :wa 1 1 
Reading  zone  data 

Erbedded  blocks  in  corrent  or  unknown:  ( 

Found  end  of  section  in  unknown:) 

.Reading  uniforr  faces 

Read  zonel:8  nare : Bottor-8  pa tchTypelD : wa 1 1 
Reading  zone  data 

Erbedded  blocks  in  corrent  or  unknown:  ( 

Found  end  of  section  in  unknown:) 

.Reading  uniforr  faces 


(screen  output  continues  on  the  next  page...) 
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Read  zonel:9  name : side_2-9  patchTypelD :wal 1 
Reading  zone  data 


FINISHED  LEX1NG 


dimension  of  grid:  3 

Creating  shapes  for  3-D  cells 

Building  patch-less  mesh... — >  FOAM  Warning  : 

From  function  polyMesh: : polyMesh( . . .  construct  from  shapes...) 
in  file  meshes/polyMesh/polyMeshFromShapeMesh.C  at  line  581 
Found  9520  undefined  faces  in  mesh;  adding  to  default  patch, 
done . 


Building  boundary  and  internal  patches. 


Creating  patch 
Creating  patch 
Creating  patch 
Creating  patch 
Creating  patch 
Creating  patch 
Creating  patch 


for  zone: 
for  zone: 
for  zone: 
for  zone: 
for  zone: 
for  zone: 
for  zone: 


3  start:  1  end:  9112  type:  interior  name:  interior-3 

4  start:  9113  end:  9180  type:  wall  name:  Inlet-4 

5  start:  9181  end:  9248  type:  wall  name:  Top-5 

6  start:  9249  end:  13872  type:  wall  name:  side_l-6 

7  start:  13873  end:  13940  type:  wall  name:  Outlet-7 

8  start:  13941  end:  14008  type:  wall  name:  Bottom-8 

9  start:  14009  end:  18632  type:  wall  name:  side_2-9 

Patch  interior-3  is  internal  to  the  mesh  and  is  not  being  added  to  the  boundary. 

Adding  new  patch  Inlet-4  of  type  wall  as  patch  0 

Adding  new  patch  Top-5  of  type  wall  as  patch  1 
Adding  new  patch  side__l-6  of  type  wall  as  patch  2 

Adding  new  patch  Outlet-7  of  type  wall  as  patch  3 

Adding  new  patch  Bottom-8  of  type  wall  as  patch  4 

Adding  new  patch  side_2-9  of  type  wall  as  patch  5 

Default  patch  type  set  to  empty 

Writing  mesh...  to  ’’constant/polyMesh"  done. 


End 


constant/  directory  and  the  createPatch  Command 

All  of  the  mesh  geometry  details  are  stored  in  the  constant/polyMesh  directory.  The 
boundary  file  is  typically  the  only  one  in  polyMesh  directory  that  gets  edited. 

Now  open  up  your  constant/polyMesh/boundary  file,  which  contains  all  of  the 
information  for  surfaces  that  were  imported  from  your  mesh.  It  should  look  like  the  screen 
capture  seen  below. 

The  information  at  the  top  of  the  file  (under  the  FoamFile  header)  gives  a  description  of 
the  file  (version,  format,  class,  etc.)  and  is  most  likely  only  useful  to  more  experienced  users,  but 
at  the  very  least  it  is  always  useful  as  a  label  to  let  you  know  where  you  are.  Further  down  the 
file  you  can  see  that  6  surfaces  (Inlet-4,  Top-5,  etc.)  were  imported  with  the  same  boundary 
names  that  were  created  in  Gridgen.  Each  surface  is  described  by  the  type  of  OpenFoam  surface, 
nFaccs  and  startFace.  Only  the  type  is  of  concern  to  the  user  at  this  stage  and  that  will  be 
discussed  in  more  detail  later  on. 


64 


w 

/ 

F  ield 

1 

|  OpenFOAM:  The  Open  Source  CFD  Toolbox 

w 

/ 

0  peration 

|  Version:  1.5-dev 

\\  / 

A  nd 

|  Web:  http://www.0penF0AM.org 

w/ 

M  anipulation 

1 

*\ 


V 

FoanFile 

{ 


V 


version 

format 

class 

location 

object 


2.0; 

ascii ; 

polyBoundaryMesh ; 
"constant/polyMesh" ; 
boundary; 


// 


Inlet-4 

{ 


type 

nFaces 

startFace 


Top-5 

{ 


type 

nFaces 

startFace 


} 

side_l-6 

{ 

type 

nFaces 

startFace 

} 

Out  let-7 

{ 

type 

nFaces 

startFace 

} 

Bottom-8 

{ 

type 

nFaces 

startFace 

} 

side_2-9 


{ 


type 

nFaces 

startFace 


wall; 

68; 

9112; 


wall; 

68; 

9180; 


wall 

4624 

9248 


wall ; 
68; 

13872; 


wall; 

68; 

13940; 


wall; 

4624; 

14008; 


// 


//  ** 


k-************************************************** 


// 
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Often  times  the  boundary  file  needs  to  be  altered  from  what  is  originally  created  during 
the  import  process.  For  this  ease  we  need  to  edit  the  type  of  surfaee  for  some  surfaces,  and  we 
would  like  to  group  certain  surfaces  together  to  avoid  redundancy  and  make  life  easier. 

First  we  will  group  some  of  the  surfaces  together  for  case  of  book-keeping.  To  group 
surfaces  together,  wc  use  the  crcatePatch  command.  For  now  let's  say  that  we  want  to  create  a 
group  of  surfaces  for  what  will  be:  the  moving  lid  (movingWall),  the  stationary  walls 
( fixedWalls ),  and  the  2-D  surfaces  on  the  front  and  back  (frontAndBack ).  The  file  that  allows  us 
to  group  surfaces  is  called  createPatchDict  and  it  is  located  in  the  system  directory  . 

If  you  open  your  system/create PatchDict  file  you  will  notice  that  it  needs  to  be  edited 
Again,  there  is  detailed  information  about  this  dictionary  file  underneath  the  FoamFilc  header, 
and  can  be  used  as  a  reference  to  the  user.  Underneath  the  FoamFilc  section  are  the 
matchToTolerance  and  pointSync  commands  which  are  not  important  right  now  and  should  be 
left  as  is.  Next,  you  will  see  a  patches  section,  which  is  where  we  will  do  our  editing  to  join 
surfaces  under  one  boundary  name  and  type . 

The  user  should  edit  their  createPatchDict  file  patchesQ  section  to  look  like  the  screen 
capture  on  the  next  page.  The  first  patch  is  essentially  renaming  the  Top-5  boundary  to 
movingWall ,  while  leaving  the  type  as  wall.  This  is  equivalent  to  simply  changing  the  name  in 
the  original  constant/polyMesh/boundary  file,  but  was  done  here  for  educational  purposes.  The 
second  patch  will  group  the  Inlet-4,  Bottom-8,  and  Outlet-7  surfaces  into  one  boundary'  named 
fixedWalls ,  and  this  new  boundary  will  remain  a  type  wall.  Lastly,  the  side  1-6  and  side_2-9 
surfaces  will  be  combined  to  frontAndBack  and  this  new'  boundary  w  ill  be  of  type  empty.  The 
empty  patch  type  is  required  for  ALL  2-D  surfaces. 

Now  we  are  ready  to  combine  the  surfaces,  but  first  it  is  generally  a  good  idea  to  copy 
our  constant  directory  before  wc  combine  our  surfaces,  so  there  is  always  a  reference.  So  for 
Linux  users,  simply: 

»  cp  -r  constant  orig ^constant 

to  copy  our  reference  constant  folder  to  orig  eonstant.  Then  in  the  case  directory  enter 
“create Patch"  in  the  eommand  line.  The  resulting  screen  dump  should  look  like  the  following 
screen  capture. 
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OpenFOAM : 

The  Open  Source  CFD  Toolbox 

\\  / 
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peration 
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Version; 

1.0 

\\  / 
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nd 

i 

Web : 

http: //www . openf oam . org 

w/ 

M 

anipulation 

i 

FoamFile 

{ 

version 
format 

root 
case 

instance 
local 

class 
object 

} 

//  ************+************************// 

//  Tolerance  used  in  matching  faces.  Absolute  tolerance  is  span  of 
//  face  times  this  factor. 
matchTolerance  IE-3; 

//  Do  a  synchronisation  of  coupled  points. 
pointSync  true; 


2.0: 
ascii ; 

"/home/penf old/mat t ijs/foam/mat  ti js2 . 1/run/ icoFoam” ; 
"cavity” ; 

"system” ; 


dictionary ; 
createPatcheDic t ; 


//  Patches  to  create. 

//  If  no  patches  does  a  coupled  point  and  face  synchronisation  anyway, 
patches 
( 


( 

name  movingWall; 
type  wall; 

constructFrom  patches; 
patches  (  Top- 5  ); 

} 

( 

name  fixedWalls; 
type  wall; 

constructFrom  patches; 

patches  (  Inlet-4  Bottom-8  Outlet-7  ); 

} 

{ 

name  f rontAndBack ; 
type  empty; 

constructFrom  patches; 
patches  (side_l-6  side_2-9  ); 

} 

); 
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[kdelaney^Retech  allCoarseStruct ]£  cp  -r  constant/  orig_constant 
[kdelaney^Retech  allCoarseStruct ] £ 

[kdelaneytfRetech  allCoarseStruct ]£ 

[kdelaney^Retech  allCoarseStruct ] £ 

[kdelaney\*Retech  allCoarseStruct ] £ 

[kdelaneyORetech  allCoarseStruct ] £  createPatch 


w 

/ 

F  ield 

|  OpenFOAM: 

The  Open  Source  CFD  Toolbox 

w 

/ 

0  peration 

|  Version: 

1 . 5 -dev 

\\  / 
w/ 

A  nd 

M  anipulation 

|  Web: 

ht  tp : //www . OpenFOAM . org 

Exec  :  createPatch 

Date  :  Mar  09  2010 

Tine  :  16:53:16 

Host  :  Retech 

PID  :  5158 

Case  :  /hone/kdelaney/NavyFoan_runs/run/f lat_plate/nov ingWallCav i ty/al ICoarseStruct 

nProcs  :  1 


//*****♦* 

Create  tine 


*  *  // 


Reading  createPatchDict . 

Using  relative  tolerance  0.001  to  natch  up  faces  and  points 
Create  polyMesh  for  tine  =  0 

Adding  new  patch  novingWall  of  type  wall  as  patch  6 
Adding  new  patch  fixedWalls  of  type  wall  as  patch  7 
Adding  new  patch  frontAndBack  of  type  enpty  as  patch  8 

Moving  faces  fron  patch  Top-5  to  patch  6 
Moving  faces  fron  patch  Inlet-4  to  patch  7 
Moving  faces  fron  patch  Botton-8  to  patch  7 
Moving  faces  fron  patch  Outlet-7  to  patch  7 
Moving  faces  fron  patch  side_l-6  to  patch  8 
Moving  faces  fron  patch  side_2-9  to  patch  8 

Doing  topology  nodification  to  order  faces. 

Synchronising  points. 

Points  changed  by  average :0  naxrO 

Renoving  patches  with  no  faces  in  then. 

Renoving  enpty  patch  Inlet-4  at  position  0 
Renoving  enpty  patch  Top-5  at  position  1 
Renoving  enpty  patch  side_l-6  at  position  2 

Renoving  enpty  patch  Outlet-7  at  position  3 

Renoving  enpty  patch  Botton-8  at  position  4 

Renoving  enpty  patch  side_2-9  at  position  5 

Renoving  patches. 

Writing  repatched  nesh  to  0.005 

End 
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You  will  notice  that  the  second  to  last  line  states  “Writing  repatehed  mesh  to  0.005”.  The 
createPatch  command  will  write  the  new  patched  surface  information  into  a  directory'  whose 
name  is  the  first  time  step  output  of  your  future  run,  which  in  this  case  is  a  0.005. 

So  now  you  need  to  get  rid  of  the  old  constant  directory  and  move  the  new  0.005 
directory  to  constant  In  Linux  this  would  be  accomplished  by: 

»  mv  constant  orig_constant 

»  mv  0.005  constant 

Now  the  transportProperties  file  needs  to  be  placed  in  the  constant  directory.  The 
transportProperties  file  must  ALWAYS  be  present  in  the  constant  directory.  In  Linux  this 
would  be  accomplished  by: 

»  mv  transportProperties  constant/ 

Now  if  you  open  the  const  ant/poly  Mesh/boundary  file  it  should  have  the  correct  number 
of  patches,  names,  and  types.  Your  new  boundary  file  should  look  like  the  screen  capture  on  the 
next  page. 
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0*- 


\* - 

FoamFile 

{ 


C++ 


-*\ 


w 
\\  / 
\\  / 
w/ 


/ 


F  ield  |  OpenFOAM:  The  Open  Source  CFD  Toolbox 

0  peratlon  |  Version:  1.5-dev 

A  nd  |  Revision:  exported 

M  anipulation  |  Web:  http://www.OpenFOAM.org 


version 
format 
class 
locat ion 
object 


2.0; 

ascii; 

polyBoundaryMesh; 
"polyMesh" ; 
boundary; 


// 


******** 

movinglv'all 

{ 

*  *  *  *  * 

type 

wall ; 

nFaces 

68; 

startFace 

} 

f lxedWal Is 
{ 

9112; 

type 

wall ; 

nFaces 

204; 

startFace 

} 

frontAndBack 

{ 

9180; 

type 

empty 

nFaces 

9248; 

startFace 

9384; 

// 


II  ***** 


************ 


*********************** 


1 1 


Now  we  have  the  geometry  imported  and  named  as  we  want  for  the  run.  A  good  next  step 
is  to  export  the  geometry  into  a  visual  package  (EnSight,  ParaView,  ete.)  and  make  sure  that  all 
surfaces  are  grouped  and  labeled  correctly.  To  export  the  geometry,  use  foaniToEnsight  for 
EnSight,  foamToVTK  for  ParaView,  and  no  additional  command  is  needed  for  ParaFoam.  So 
now  take  a  minute  or  two  and  inspect  your  geometry  in  your  package  of  choice.  Your  geometry 
should  look  like  the  below  Figure,  with  the  appropriate  surface  labels. 


70 


Cavity  geometry  as  seen  in  EnSight 


Material  Properties 

The  next  step  is  to  set  up  the  material  properties  for  the  fluid.  For  the  icoFoam  solver, 
only  the  kinematic  viscosity  is  required  in  the  constant/transportProperties  dictionary'  file.  Open 
the  transportProperties  file,  it  should  look  like  the  screen  capture  below.  No  editing  is 
necessary,  just  note  that  kinematic  viscosity  is  always  set  in  transportProperties. 


0 


*-  C++ 


4 


I  \\  / 

I  \\  / 

I  \\  / 

I  w/ 

\* - 

FoamFile 

{ 

version 

format 

class 

object 

} 

II  *  *  *  *  + 


F  ield 
0  peration 
A  nd 

M  anipulation 


OpenFOAM:  The  Open  Source  CFD  Toolbox 
Version:  1.5 

Web:  http://www.OpenFOAM.org 


2.0; 
ascii ; 
dictionary; 
transportProperties; 

*4  +  +  +  ***  +  *****  +  ******  +  ****** 


4  4 


*\ 

I 

I 

I 


V 


*  // 


nu 


nu  [0  2  -1  0  0  0  0]  0.01; 


// 


4444444444444444444444444444444444444444444444444444444444444444444444*4* 


// 
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The  kinematic  viscosity,  nu,  is  entered  in  as: 
nu  nu  [0  2-1  0  0  0]  0.01; 

where  [02-1  000  0]  sets  the  units  based  on  [Mass  Length  Time  Temperature  Quantity  Current 
Luminous  intensity].  The  kinematic  viscosity  dimensions  are  Length /Time  or  in  SI  units:  ms/s. 
The  value  of  the  kinematic  viscosity  is  set  to  0.01.  Remember  that  this  value  must  always  be 
consistent  with  the  Reynolds  number, 

Re  =  UL/nu 

0/  directory  (Initial  and  Boundary  Conditions) 

Now  we  turn  our  attention  to  the  initial  and  boundary  conditions,  which  are  stored  in  the 
0/directory. 

For  the  icoFoam  solver  only  U  and  p  files  are  needed  in  the  0/ directory. 

Open  the  0/U  file:  it  should  look  like  the  screen  capture  below. 
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-  C++  -* 


\\  /  F  ield 

\\  /  0  peration 

\\  /  A  nd 

\\/  M  anipulation 


OpenFOAM:  The  Open  Source  CFD  Toolbox 
Version:  1.5 

Web :  ht tp : //www . OpenFOAM . org 


V - 

FoamFile 

{ 

version 

format 

class 


} 

// 


object 


2.0; 

ascii; 

volVectorField; 

U; 


*********  ************************** 


*  It 


dimensions  [01-10000]; 

internalField  uniform  (000); 


boundaryField 

{ 

movingWall 

{ 

type 

value 

} 


f ixedWalls 

{ 

type 

value 

} 

frontAndBack 

{ 

type 

} 


f ixedValue ; 
uniform  (100); 


f ixedValue ; 
uniform  (000): 


empty; 


//  +  ************************  ***♦*****************************************,; 


// 


In  the  0/U  file  notice: 

U  is  a  vector  field  quantity.  All  U  values  must  be  set  in  vector  format,  (XX X). 

The  U dimensions  must  match  the  variable  by  so  velocity  is  [0  1  -/  0  0  0  0] 

The  ‘fintemalField,’  sets  the  initial  flow  field  condition  for  U.  For  this  case  the  fluid  is 
initially  at  rest  (Ux,Uy,Uz  =  0) 
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The  “boundary Field'’  sets  velocity  boundary  conditions  for  ALL  surfaces.  All  surfaces 
must  be  included  with  proper  BC’s.  All  surface  names  must  match  the 
constant/polyMesh/boundary  surface  name  exactly.  The  order  of  surfaces  is  not  important,  but 
the  names  must  match  identically. 

The  velocity  boundary  conditions  for  the  three  surfaces  are  as  follows: 

The  movingWall  is  set  with 

type  fixedValue; 
value  uniform  (1  0  0); 
for  the  top  lid  to  move  with  velocity  of  Ux=\. 

The  fixed  Walls  are  set  with 

type  fixedValue; 
value  uniform  (0  0  0); 

to  apply  the  no-slip  boundary  condition  to  the  walls. 

The  frontAndBack  surfaces  are  set  with 
type  empty; 

because  they  are  the  2-D  boundaries.  ALL  2-D  boundaries  must  have  type  set  to  empty. 

Again,  the  order  of  the  surfaces  in  the  0A..  files  doesn’t  matter,  but  the  names  and  types 
MUST  be  consistent  with  those  listed  in  the  constant/polyMesh/boundary >  file. 

Now  open  the  0/p  file,  it  should  look  like  the  screen  capture  below. 
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-  C++  - 


w 

/ 

F  ield 

1 

|  OpenFOAM: 

The  Open  Source  CFD  Toolbox 

w 

/ 

0  peration 

|  Version: 

1.5 

\\  / 

A  nd 

|  Web: 

ht  tp :  //www .  OpenFOAM .  or  g 

w/ 

M  anipulation 

1 

v - V 

FoamFile 

{ 

version  2.0; 

format  ascii; 

class  volScalarField; 

object  p; 

} 

yy  *************************************// 

dimensions  [02-20000]; 

internalField  uniform  0; 


boundaryField 

{ 

movingV/all 

{ 

type  zeroGradient ; 

} 

fixedWalls 

{ 

type  zeroGradient; 

} 

frontAndBack 

{ 

type  empty; 

} 


} 


// 


// 


In  the  0/p  file  notice: 

p  is  a  scalar  field  quantity.  All  p  values  must  be  set  as  a  sealar,  X ,  value. 

The  p  dimensions  must  mateh  the  variable  by  M,L,T,...  so  pressure  is  [0  2  -2  0  0  0  0], 
because  in  icoFoam  p  is  actually  the  pressure  divided  by  the  density,  thus  the  SI  units  would  be 

itT/s*. 

The  “internalField'’  sets  the  initial  condition  for  p.  For  this  case  (and  most  others)  we  do 
not  care  about  the  absolute  value  of  the  pressure,/?,  so  we  just  set  it  to  0  for  ease. 
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The  boundaryFicld  sets  pressure  boundary  conditions  for  ALL  surfaces.  All  surfaces 
must  be  included  with  proper  BC’s  that  are  consistent  with  the  constant/polyMcsh/boundary 
surfaces. 

The  pressure  boundary  conditions  for  the  three  surfaces  are  as  follows: 

The  movingWall  is  set  with 

type  zeroGradient; 

for  the  moving  wall. 

The  jixedWalls  arc  set  with 

type  zeroGradient; 

for  the  no-slip  wall. 

The  frontAndBack  surfaces  are  set  with 
type  empty; 

because  they  are  the  2-D  boundaries.  ALL  2-D  boundaries  must  have  type  set  to  empty, 
system/  directory  (Solver  Settings) 

Now  we  will  look  at  some  of  the  solver  settings  and  controls  that  are  located  in  the 
system/  directory.  We  will  focus  on  the  eontrolDict ,  fvSolution,  and  fvSchemcs  files.  We 
already  used  the  createPatchDict  to  merge  multiple  surfaces. 

Open  the  system/ eontrolDict  dictionary  file.  It  should  look  like  the  screen  capture  below. 
The  eontrolDict  file  sets  all  of  the  run-time  parameters  and  output  information.  This  is  also 
where  run-time  libraries  and  functions,  such  as  force  outputs  over  a  patch  and  dynamic  mesh 
libraries  are  specified. 
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0' - *-  C++  -* - -  \ 


w 

/  F  ield 

|  OpenFOAM:  The  Open  Source  CFD  Toolbox 

w 

/  0  peration 

|  Version;  1.5 

\\  / 

A  nd 

|  Web:  http://www.OpenFOAM.org 

w/ 

M  anipulation 

FoarFi  le 


{ 

version 

fornat 

class 

object 

} 

//***** 


2.0; 
ascii ; 
dictionary; 
controlDict ; 


*********  ii 


application  icoFoati; 


startFron 

startTine ; 

startTine 

0; 

stopAt 

endTine ; 

endTine 

0.3; 

deltaT 

0.00001; 

writeControl 

t ineStep ; 

writelnterval 

3000; 

purgeWrite 

0; 

wri teFornat 

ascii ; 

writePrecision 

6; 

writeConpression  uncompressed; 
tineFornat  general; 

t inePrecision  6; 

runTineModif iable  yes; 


//  *************************************************************************  i f 


I- 


The  solver  specified  in  application  input  does  not  matter.  The  solver  is  specified  on  the 
command  line  or  in  a  script  file.  Thus,  this  is  an  insignificant  line  for  our  purposes. 

The  solver  settings  are  fairly  obvious,  and  more  detail  is  provided  on  page  U-108  of  the 
User’s  Guide  (http://foam.sourceforge.net/doc/Guides-a4/UserGuidc.pdf).  For  now  we  will  only 
cover  a  broad  view  of  the  file. 
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Wc  know  that  icoFoam  is  a  transient  solver.  We  see  that  we  will  run  the  simulation  from 
0  to  0.3  seconds  in  time  steps  of  0.00001  seconds  with  the  solver  writing  output  information 
every  3000  time  steps  (t=0.03,  0.06,  0.09,  ...).  The  data  will  be  written  in  ASCII  format  in 
directories  that  are  denoted  by  time  6  digits  long.  Notice  that  runTimeModifiable  is  chosen  to 
yes ,  this  means  that  we  can  make  changes  to  the  controlDict  in  the  middle  of  a  run,  and  they  will 
be  adjusted  on  the  fly,  as  opposed  to  having  the  settings  set  in  stone  for  the  whole  calculation. 

One  important  note  is  that  to  start  a  calculation  from  a  previous  solution  the  stariFrom 
entry  must  be  switched  to  latestTime ,  and  desired  start  time  information  (directory  and  BC's) 
must  be  present  in  the  ease  directory. 

Now  open  the  system/fvSolution  dictionary  file.  It  should  look  like  the  screen  capture  on 
the  next  page. 

The  fvSolution  file  contains  linear  solver  information  as  well  as  solver  algorithm  settings. 

The  solvers  section  contains  linear  solver  settings  for  pressure  and  velocity.  Note  that  for 
this  ease  we  arc  using  preconditioned  conjugate  gradient  solvers  (PCG  for  symmetric  matrices 
and  PBiCG  for  asymmetric  matrices),  but  we  also  commonly  use  multi-grid  solvers  (GAMG, 
AAMG,  etc.).  The  solver  tolerance  and  relative  tolerance  settings  are  not  important  right  now. 
The  minlter  command  sets  a  minimum  number  of  times  the  linear  solver  will  iterate  on  a 
variable.  It  is  usually  recommended  that  the  user  always  set  a  minimum  number  of  iterations  >  0 
to  prevent  the  solver  from  prematurely  not  solving  for  a  variable. 

Below  the  solvers  section  arc  pressure-implicit  split-operator  (PISO)  algorithm  control 
settings.  These  PISO  settings  are  not  particularly  useful  to  the  user  at  this  time,  so  only  a  broad 
view  of  what  each  setting  means  is  given.  Also,  note  that  the  PISO  algorithm  must  be  used  for 
all  transient  solvers  and  the  SIMPLE  algorithm  must  be  used  for  all  steady-state  solvers.  For  this 
case  wc  have  nCorrectors  set  to  2,  which  means  that  we  will  solve  the  pressure  equation  twice 
per  time  iteration.  The  value  of  nNonOrthogonalCorrectors  is  set  to  0.  This  parameter  is  not 
particularly  important  to  the  user  at  this  moment.  Notice  that  wc  have  set  cell  number  0  as  our 
reference  cell,  where  the  reference  value  is  0.  This  is  the  reference  pressure  for  the 
incompressible  solver. 
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0*- 


-  C++ 


\\  /  F  ield 

\\  /  0  peration 

\\  /  A  nd 

\\/  M  anipulation 


OpenFOAM:  The  Open  Source  CFD  Toolbox 
Version:  1.5 

Web :  http://www.OpenFOAM.org 


\* - 

FoamFile 

{ 

version 

format 

class 

object 

} 

//****+ 


2.0; 
ascii ; 
dictionary ; 
fvSolution; 


*  *  // 


solvers 

{ 

p  PCC 

{ 

preconditioner 

tolerance 

relTol 

minlter 

}; 


DIC; 
le-06 ; 
0; 

5; 


U  PBiCG 

{ 

preconditioner  DILU; 
tolerance  le-05; 

relTol  0; 

minlter  5; 

}; 


PISO 

{ 


nCorrectors 

2; 

nNonOrthogona 

ICorrectors  0; 

pRefCell 

0; 

pRef Value 

0; 

// 


*+***********-*************<********-*******************-*♦**-*♦***-******♦**** 


// 


Now  open  the  system/fvSchemes  dictionary  file.  It  should  look  like  the  screen  capture 

below. 

Many  of  the  fvSchemes  settings  are  not  particularly  useful  to  the  user  at  this  time,  so  only 
a  broad  view  of  what  each  setting  means  is  given.  For  more  detail  on  these  settings  consult  page 
U-l  10  of  the  User’s  Guide.  More  detail  of  the  discretization  settings  is  given  in  the  sintpleFoam 
and  rasInterFoam  tutorials. 
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*  * 

\\  / 

F  ield 

i 

I  OpenFOAM: 

The  Open  Source  CFD  Toolbox 

\\  / 

0  peration 

|  Version: 

1.5 

\\  / 

A  nd 

I  Web: 

http : //www . OpenFOAM . org 

w 

M  anipulation 

1 

V 

FoanFile 

{ 

version 

format 

class 

object 

} 

//****** 

ddtSchenes 

{ 

default 

} 

gradSchenes 

( 

default 

grad(p) 

} 

divSchemes 

{ 

default 
div(phi , U) 

} 


2.0; 
ascii ; 
dictionary ; 
fvSchemes; 


******* 


// 


Euler ; 


Gauss  linear; 
Gauss  linear; 


none; 

Gauss  linear; 


laplacianSchemes 

{ 

default  none; 

laplacian(nu,U)  Gauss  linear  corrected; 
laplacian( (1| A(U) ) ,p)  Gauss  linear  corrected; 

} 

interpolat ionSchemes 

{ 

default  linear; 

interpolate(HbyA)  linear; 

} 


snGradSchemes 

{ 

default 

} 

f luxRequired 

{ 

default 

Pi 

} 


to 


**************4 


corrected; 


no ; 


t*A********4r 


£r********************4 


************ 


*******  II 
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The  fvSchemes  file  sections  declares  the  following  settings: 
ddt  ->  time  discretization 
gradSchemes  4  gradient  term  discretizations 
divSchemes  divergence  terms  discretization 

laplacianSchemes  Laplaeian  terms  discretization 

interpolationSchemcs  interpolation  of  values  from  cell  centers  to  cell  faee  eenters 

snGradSchemes  surface  normal  gradient  evaluation  at  eell  faces 

JliixRequircd  4  lists  fields  where  flux  is  generated  in  the  application 

For  all  of  the  fvSchemes  fields  a  default  value  can  be  specified,  or  default  can  be  set  to 
none  which  means  that  the  user  must  enter  all  values  for  the  appropriate  variables  themselves. 

Running  icoFoam 

Now  we  are  ready  to  run.  Type  icoFoam  on  the  command  line  like  the  screen  capture 
below,  and  hit  4tc/r/+c”  to  take  a  look  at  the  first  few'  iterations  (piping  the  screen  dump  to  a  log 
file  by  using  “ icoFoam  \  tee  log  ’  is  another  option). 
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[kdelaneyORetech  tut_allCoarseStruct  ]  t*  icoFoan 


\\  / 

F  ield 

1 

I  OpenFOAM: 

The  Open  Source  CFD  Toolbox 

\\  / 

0  peration 

|  Version: 

1 . 5-dev 

\\  / 

A  nd 

1  Web: 

http:/ / www . OpenFOAM . org 

w/ 

H  anipulation 

1 

Exec  :  icoFoan 

Date  :  Jun  30  2010 

Tire  :  14:53:06 

Host  :  Retech 

P1D  :  20508 


Case  :  /hone/kdelaney/NavyFoaii_runs/run/f lat_plate/iK>vingWallCavity/tut_allCoarseStruct 
nProcs  :  1 


„•••••***■••••****•*•***•••*****•*•***// 
Create  tire 


Create  resh  for  tire  r  o 


Reading  transport  Proper  ties 
Reading  field  p 
Reading  field  U 

Reading/calculat ing  face  flux  field  phi 


Starting  tire  loop 
Tire  -  le-OS 

Courant  Nurber  rean:  0  rax:  0  velocity  magnitude-  0 

PBiCG:  Solving  for  Ux,  Initial  residual  =•  1,  Final  residual  =•  1.68372e-20,  No  Iterations  5 
PBiCC:  Solving  for  Uz:  solution  singularity 

PCC:  Solving  for  p,  Initial  residual  =  1,  Final  residual  =  6.77855e-07.  No  Iterations  138 
tire  step  continuity  errors  :  sur  local  =  1.14833e-15.  global  =  3.92084e-26,  curulative  =  3.92084e-26 
PCC:  Solving  for  p.  Initial  residual  =  0.0407237,  Final  residual  =  8.83925e-07,  No  Iterations  126 
tire  step  continuity  errors  :  sur  local  =  3.57796e-15,  global  =  -5.09959e-25,  curulative  =  -4.7075e-25 
Execut lonTire  -  0.13  s  ClockTire  -0s 

Tire  -  2e-0S 

Courant  Nurber  rean:  1.36506e-07  rax:  9.7988e-05  velocity  ragnitude:  0.0362317 

PBiCC:  Solving  for  Ux,  Initial  residual  =  0.134603,  Final  residual  =  1.07427e-21,  No  Iterations  5 

PBiCC:  Solving  for  Uz,  Initial  residual  =  0.333103,  Final  residual  =  8.02106e-21.  No  Iterations  5 

PCC:  Solving  for  p,  Initial  residual  =  0.100492,  Final  residual  =  8.49585e-07,  No  Iterations  125 
tire  step  continuity  errors  :  sur  local  =  3.78654e-15,  global  =  3.849e-25.  curulative  =  -8.S8S04e  26 
PCC:  Solving  for  p.  Initial  residual  =■  0.0020439,  Final  residual  =  9.96011e-07.  No  Iterations  112 
tire  step  continuity  errors  :  sur  local  =  4.8709e-15,  global  =  4.60384e-25.  curulative  =  3.74534e-25 
Execut ionTire  =  0.2  s  ClockTire  =  0  s 

Tire  =  3e-0S 

Courant  Nurber  rean:  2.66818e-07  rax:  0.000190662  velocity  ragnitude:  0.0705446 

PBiCC:  Solving  for  Ux,  Initial  residual  =  0.0711115,  Final  residual  =  5.11184e-22,  No  Iterations  5 

PBiCC:  Solvinc  for  Uz.  Initial  residual  =•  0.196949.  Final  residual  =  5.20581e-21.  No  Iterations  5 

Some  observations  from  the  first  few  iterations: 

You  can  see  that  the  solver  started  from  time  equal  to  0  seconds  and  is  marching  in 
increments  of  le-5  seconds. 

For  eaeh  iteration  the  pressure  equation  is  solved  twice  and  the  velocity  equations  are 
solved  once.  For  each  variable  linear  solver  we  can  see  the  initial  residual,  final  residual,  and  the 
number  of  iterations  it  took  to  drop  from  the  initial  to  the  final  residual.  We  set  all  of  these 
toleranees  and  iteration  criteria  in  the  system/fvSolution  dietionary  file. 

There  are  also  Courant  number  and  continuity  error  reports. 
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The  best  way  to  typically  monitor  the  solution  is  to  make  sure  that  the  velocity  magnitude 
stays  at  a  reasonable  number,  and  make  sure  that  initial  pressure  residuals  are  decreasing  or  are 
holding  steady  at  an  acceptable  value. 

The  last  line  of  the  time  iteration  produces  execution  and  clock  time  information.  This  is 
useful  in  gauging  the  efficiency  of  your  solution. 

Now  let  the  icoFoam  solver  go  until  0.5  seconds  to  get  a  converged  solution. 

Post-Processing 

Notice  that  there  are  many  time  directories  in  your  ease  directory.  Each  of  these 
directories  contains  output  information  for  their  respective  time  step. 

To  look  at  the  post-processed  results  simply  type  the  following  commands,  depending  on 
the  post-processing  tool  of  choice: 

» foamToEnSight  - IatestTime  to  look  at  the  results  in  EnSight 

>>  foamToVTK  - IatestTime  ->  to  look  at  the  results  in  ParaView 

where  the  command  ~ IatestTime  is  used  to  only  look  at  the  results  from  the  last  output  time  step. 
To  look  at  the  results  for  all  time  steps  simply  leave  off  the  -IatestTime  command,  and  to  look  at 
the  results  for  a  specific  time  (ie  0.005)  us  e-time  0.005. 

To  look  at  the  results  in  ParaFoam,  no  additional  commands  are  needed,  simply  open 
ParaFoam  in  the  case  directory. 

Your  results  should  look  like  the  velocity  magnitude  (Vmag)  and  axial  velocity 
( VmagfxJ )  contours  below. 
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Appendix  D:  simpleFOAM  Body-1  Tutorial 

This  tutorial  involves  using  the  turbulent,  steady,  incompressible  solver  for  a  3-D  body- 
of-revolution,  the  Body-1.  Only  half  the  body  is  solved,  as  symmetry  is  assumed.  The  domain  is 
non-dimensionalized  by  length,  so  all  lengths  in  the  domain  are  normalized  by  body  length.  First, 
we  will  go  over  pre-processing  and  case  setup,  then  we  will  run  the  test  case,  and  Finally  we  will 
look  at  some  post-processed  results. 

For  more  detailed  information  on  the  OpenFOAM  code  and  settings  consult  the  User's 
Guide:  http://foam.sourceforge.net/doc/Guides-a4/UserGuide.pdf. 

Pre-Processing  and  Case  Setup 

Upon  looking  in  the  ease  directory  (sereen- eapture  below)  you  will  notice  that  there  are 
directories  labeled  system  and  0 ,  and  a  file  named  transportPropcrties  There  is  a  parallel 
processing  script  named  oFOAM.scp .  There  is  also  a  file  labeled  bodyl  Box_ASCILfluentcas , 
which  is  a  mesh  exported  from  Gridgen  in  Fluent  ASCII  double  precision  format.  We  will 
discuss  the  existing  Files  and  directories  later,  but  now  we  need  to  import  the  fluent  .cas  File  into 
OpenFoam.  This  will  be  done  using  the  command  JluentMeshToFoam . 

[delaneykflamazon  bodyl_Tutorial] £  1 
total  237M 


-rw-r — r — 

1 

delaneyk 

users 

1.6K 

Apr 

8 

15:17 

transportPropert ies 

-rwxr-xr-x 

1 

delaneyk 

users 

4.6K 

Apr 

8 

15:17 

RASProperties 

drwxr-xr-x 

2 

delaneyk 

users 

80 

Apr 

8 

15:17 

system 

drwxr-xr — 

2 

delaneyk 

users 

81 

Apr 

8 

15:17 

0 

-rw-r — r — 

1 

delaneyk 

users 

658 

Apr 

8 

15:17 

oFOAM. sep 

-rw-r--r-- 

1 

delaneyk 

users 

237M 

Apr 

8 

15:31 

bodyl_Box-ASCI I . f luent .cas 

Mesh  Input 

Now  enter  the  JluentMeshToFoam  command  as  seen  below  and  view  the  output  from 
the  screen  dump.  Your  output  should  be  the  same  as  the  screen  captures  on  the  next  pages. 

There  is  a  lot  of  information  on  the  screen  dump,  most  of  which  is  self  explanatory.  The 
most  important  part  to  notice  is  the  last  two  lines  of  text  in  the  screen  dump  which  tell  you  that 
the  mesh  information  has  been  written  into  a  directory  named  polyMesh  inside  a  newly  created 
directory  named  constant .  Finally  the  command  ends  successfully  with  the  “End.”  statement 
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[delaneyk^anazon  bodyl_Tutorial ]$  f IuentMeshToFoan  body l_Box-ASCII . fluent .cas 

/* 


\\  / 

F  ield 

1 

1 

OpenFOAM : 

The  Open  Source  CFD  Toolbox 

\\  / 

0  peration 

1 

Version: 

1 . 5-dev 

\\  / 

A  nd 

1 

Revision: 

exported 

w/ 

M  anipulation 

1 

Web: 

ht  tp : //www . OpenFOAM . org 

V 


\* 

Exec  :  f luentWeshToFoan  body l_Box-ASCII . fluent . cas 

Date  :  Apr  08  2010 

Tine  :  15:31:37 

Host  :  anazon.dt .navy. nil 

PID  :  19865 

Case  :  /san/hone/delaneyk/NavyFOAM-1 . 5-dev-rev995/delaneyk-l . 5-dev/run/bodyl/boundingBoxFar 

nProcs  :  1 


//•*•** 

Create  tine 


*  *  *  // 


-->  FOAM  Warning  : 

Fron  function  dlLibraryTable: :open(const  fileNaneA  funct ionLibNane) 
in  file  db/dILibraryTable/dlLibraryTable.C  at  line  86 

could  not  load  /san/hone/de laneyk’/NavyFOAM-l . 5-dev-rev995/Navy FOAM/1 ib/linux64CccDPOpt  H i 
Model llpr intCoef fsEv 
Dinension  of  grid:  3 
■MNunber  of  cells:  2180547 
nunber  of  faces:  5106844 
Nunber  of  points:  860613 
Other  readCellCroupData :  a  1  2145c3  1  0 
Reading  nixed  cells 
Reading  nixed  faces 
Reading  uniforn  faces 
Reading  uniforn  faces 
Reading  uniforn  faces 
Reading  uniforn  faces 
Reading  uniforn  faces 
Reading  uniforn  faces 
Reading  nixed  faces 
Reading  points 

Read  zone2:10  nane: fluid  patchTypelD : fluid 
Reading  zone  data 

Read  zone2:8  nane : synne try  patchTypelD : synne try 
Reading  zone  data 

Read  zone2:6  nane: out  let  patchTypelD :pressure-out let 
Reading  zone  data 

Read  zone2:5  nane: inlet  patchTypelD: inlet-vent 
Reading  zone  data 

I 

Read  zone2:4  nane : f arf ield  patchTypelD :pressure-f ar-f ield 
Reading  zone  data 

Read  zone2:l  nane: bow  patchTypeID:wal I 
Reading  zone  data 

Read  zone2:2  nane:nidbody  patchTypelD : wall 
Reading  zone  data 

Read  zone2:3  nane: stern  DatchTvoelDiwall 

(screen  output  continues  on  the  next  page... ) 
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Reading  zone  data 

Read  zone2:9  name : interior-faces  patchTypelD: interior 
Reading  zone  data 


FINISHED  LEXING 


dimension  of  grid:  3 

Creating  shapes  for  3-D  cells 

Building  patch-less  mesh... — >  FOAM  Warning  : 

From  function  polyMesh: :polyMesh( . . .  construct  from  shapes...) 
in  file  meshes/polyMesh/polyMeshFromShapeMesh.C  at  line  581 
Found  111916  undefined  faces  in  mesh;  adding  to  default  patch, 
done . 


Building  boundary  and  internal  patches. 
maxZonelD:  9 

for  zone 
for  zone 
for  zone 
for  zone 
for  zone 
for  zone 
for  zone 
for  zone 


Creating  patch 
Creating  patch 
Creating  patch 
Creating  patch 
Creating  patch 
Creating  patch 
Creating  patch 
Creating  patch 


1  end:  4994928  type: 
4994929  end:  5001640 
5001641  end: 

5013073  end: 

5035505  end: 

5037283  end: 

5037725  end: 

5038167  end: 


5013072 

5035504 

5037282 

5037724 

5038166 

5106844 


start : 
start : 
start : 
start : 
start : 
start : 
start : 

8  start: 

Patch  interior-faces  is  internal  to  the  mesh  and  is  not 
Adding  new  patch  bow  of  type  wall  as  patch  0 
Adding  new  patch  midbody  of  type  wall  as  patch  1 
Adding  new  patch  stern  of  type  wall  as  patch  2 
Adding  new  patch  farfield  of  type  patch  as  patch  3 
Adding  new  patch  inlet  of  type  patch  as  patch  4 
Adding  new  patch  outlet  of  type  patch  as  patch  5 
Adding  new  patch  symmetry  of  type  symnetry Plane  as  patch 


interior  name:  interior-faces 
type:  wall  name:  bow 
type:  wall  name:  midbody 
type:  wall  name:  stern 
type:  pressure-far-f ield  name:  farfi 
type:  inlet-vent  name:  inlet 
type:  pressure-outlet  nane:  outlet 
type:  symmetry  name:  symmetry 
being  added  to  the  boundary. 


Default  patch  type  set  to  empty 

Writing  mesh...  to  "const ant/polyMeshM  done. 


End 


Now  run  the  checkMesh  command  for  two  reasons: 

•  to  make  sure  the  mesh  was  imported  eorreetly 

•  to  asses  the  quality  of  the  mesh  for  the  OpenFOAM  solver 

Your  checkMesh  output  should  look  like  the  screen  captures  on  the  next  pages. 
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[delaneyk£amazon  bodyl_Tutorial]£  checkMesh 


/ 


\ 


w  / 

F 

ield 

1 

|  OpenFOAM:  The  Open  Source  CFD  Toolbox 

\\  / 

0 

peration 

|  Version:  1.5-dev 

\\  / 

A 

nd 

|  Revision:  exported 

\\! 

M 

anipulation 

|  Web:  http://www.0penF0AM.org 

Exec  :  checkMesh 

Date  :  Apr  08  2010 

Tine  :  15:36:20 

Host  :  amazon.dt .navy .nil 

PID  :  21145 


*\ 


V 


Case  :  /san/hone/delaneyk/NavyFOAM-1 . 5-dev-rev995/delaneyk-l . 5-dev/run/bodyl/boundingBo\Far 
nProcs  :  1 


//*************♦♦**  *  ****************  *  *  *  // 
Create  tine 


— >  FOAM  Warning  : 

Fron  function  dlLibraryTable : : open(const  fileNaneA  funct ionLibNane) 
in  file  db/dlLibraryTable/dlLibraryTable.C  at  line  86 

could  not  load  /san/hone/delaneyk/NavyFOAM-1 . S-dev-rev995/NavvFOAM/lib/linux64CccDPOpt /li 
Model llprintCoef f sEv 
Create  polyMesh  for  tine  =  constant 

Tine  =  constant 


Mesh  stats 

points:  860613 

faces:  5106844 

internal  faces:  4994928 

cells:  2180547 

boundary  patches:  7 

point  zones:  0 

face  zones:  0 

cell  zones:  0 


Nunber  of  cells 
hexahedra : 
pnsns : 
wedges : 
py ran ids: 
tet  wedges: 
tetrahedra : 
polyhedra : 


of  each  type: 
0 

1379584 

0 

0 

0 

800963 

0 


Checking  topology. . . 

Boundary  definition  OK. 

Point  usage  OK. 

Upper  triangular  ordering  OK. 
Face  vertices  OK. 

Number  of  regions:  1  (OK). 


Checking  patch  topology  for  nultiply  connected  surfaces  . . . 


Patch 

Faces 

Points 

Surface  topology 

bow 

6712 

3459 

ok 

(non-closed  singly  connected) 

nidbody 

11432 

5882 

ok 

(non-closed  singly  connected) 

stern 

22432 

11490 

ok 

(non-closed  singly  connected) 

farf ield 

1778 

950 

ok 

(non-closed  singly  connected) 

inlet 

442 

252 

ok 

(non-closed  singly  connected) 

outlet 

442 

252 

ok 

(non-closed  singly  connected) 

(Screen  capture  continues  on  next  page...  ) 
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symmetry 


68678 


50479 


ok  (non-closed  singly  connected) 


Checking  geometry... 

This  is  a  3-D  mesh 

Overall  domain  bounding  box  (-10  -10  0)  (10  10  10) 

Mesh  (non-empty)  directions  (111) 

Mesh  (non-empty,  non-wedge)  dimensions  3 

Boundary  openness  (-1 . 71587e-18  -1.2l515e-17  -5.8041e-16)  Threshold  =  le-06  OK. 

Max  cell  openness  =  5.13l85e-14  OK. 

Max  aspect  ratio  =  564.729  OK. 

Minumum  face  area  =  S.88008e-10.  Maximum  face  area  =  1.24594.  Face  area  magnitudes  OK. 
Min  volume  =  6.96628e-14.  Max  volume  =  0.421459.  Total  volume  =  4000.  Cell  volumes  OK. 
Mesh  non-orthogonality  Max:  61.1137  average:  9.70308  Threshold  =  70 
Non-orthogonality  check  OK. 

Face  pyramids  OK. 

Max  skewness  =  1.25241  OK. 

Mesh  OK. 


There  is  a  lot  (probably  too  much)  of  information  with  the  checkMesh  screen  dump.  At 
the  top,  the  Mesh  stats  section  shows  that  the  mesh  has  2.1 8  million  total  cclls/clemcnts,  and  the 
Number  of  cells  of  each  type  section  shows  1.38  million  arc  prisms  and  0.8  million  are 
tetrahedral  elements.  Below  that,  we  see  that  the  topology  checks  out  OK  and  that  all  the  surfaces 
are  correctly  connected. 

Finally,  the  Checking  geometry ...  section  displays  mesh  quality  statistics.  This  section 
gives  a  lot  of  information,  but  the  most  important  parts  are  the  aspect  ratio,  non-orthogonality, 
and  skewness.  For  this  case  all  check  out  OK ,  so  we  are  free  to  proceed  knowing  the  mesh  is  of 
high  quality. 

Sometimes  it  is  not  possible  to  create  a  mesh  without  any  high  aspect  ratio,  non- 
orthogonal,  or  skewed  cells.  In  fact,  most  meshes  created  will  contain  bad  cells,  and  run  fine. 
However,  at  some  point  (which  is  not  quantitatively  clear)  the  mesh  will  be  so  poor  it  either 
won’t  run,  or  it  will  take  a  long  time  to  run.  There  aren’t  exact  guidelines  on  OpcnFOAM  mesh 
quality;  it  simply  takes  experience  running  various  meshes. 

constant/  directory  and  the  createPatch  Command 

All  of  the  mesh  geometry'  details  arc  stored  in  the  constant/poly  Mesh  directory.  The 
boundary  file  is  typically  the  only  one  in  polyMesh  directory  that  gets  edited. 

Now  open  up  your  constant/polyMesh/boundary  file,  which  contains  all  of  the 
information  for  surfaces  that  were  imported  from  your  mesh.  It  should  look  like  the  screen 
capture  seen  below. 

The  information  at  the  top  of  the  file  (under  the  FoamFile  header)  gives  a  description  of 
the  file  (version,  format,  class,  etc.)  and  is  most  likely  only  useful  to  more  experienced  users,  but 
at  the  very  least  it  is  always  useful  as  a  label  to  let  you  know  where  you  are.  Further  down  the 
file  you  can  see  that  7  surfaces  (bow,  midbody,  etc.)  were  imported  with  the  same  boundary 
names  that  were  created  in  Gridgen.  Each  surface  is  described  by  the  type  of  GpcnFoam  surface, 
nFaces  and  startFace.  Only  the  type  is  of  concern  to  the  user  at  this  stage  and  that  w  ill  be 
discussed  in  more  detail  later  on. 
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*-  C++  - 


/ 


I  w 

i  w  / 
i  w  / 

I  w/ 

V - 

FoamFile 

{ 

version 
format 
c  lass 
location 
object 

) 

//*♦*** 


F  ield  |  OpenFOAM:  The  Open  Source  CFD  Toolbox 

0  per at ion  |  Version:  1.5-dev 

A  nd  |  Revision:  exported 

M  anipulation  |  Web:  http://www.OpenFOAM.org 


2.0; 

ascii; 

polyBoundaryMesh ; 

"constant/polyMesh" 

boundary; 


bow 

{ 


type 

nFaces 

startFace 

) 

midbody 

{ 

type 

wall ; 

6712; 

4994928; 

wall ; 

11432; 

5001640; 

nFaces 

startFace 

) 

stem 

{ 

type 

wall; 

22432; 

5013072; 

nFaces 

startFace 

) 

farf ield 
{ 

type 

patch ; 

nFaces 

1778; 

startFace 

) 

inlet 

{ 

5035504; 

type 

patch; 

nFaces 

442; 

startFace 

} 

outlet 

{ 

5037282; 

type 

patch; 

nFaces 

442; 

startFace 

) 

symmetry 

f 

5037724; 

l 

type 

symmetryPlane 

nFaces 

68678; 

startFace 

5038166; 

// 


C 

Often  times  the  boundary  file  needs  to  be  altered  from  what  is  originally  created  during 
the  import  process.  For  example,  the  hull  might  be  imported  as  5  different  surfaces  and  you 


90 


would  like  to  group  them  together  as  one  surface.  For  this  case  we  will  group  the  separate  hull 
surfaces  together  to  avoid  redundancy  and  make  life  easier. 

To  group  separate  surfaces  together,  we  use  the  crcatePatch  command.  For  now  let’s  say 
that  we  want  to  create  a  group  of  surfaces  for  what  will  be  the  hull.  The  file  that  allows  us  to 
group  surfaces  is  called  createPatchDict  and  it  is  located  in  the  system  directory. 

If  you  open  your  system/createPatchDict  file  you  will  notice  that  it  needs  to  be  edited  to 
group  the  surfaces  from  constant/poly  Mesh/boundary.  Again,  there  is  detailed  information  about 
this  dictionary  file  underneath  the  FoamFile  header,  and  can  be  used  as  a  reference  to  the  user 
Underneath  the  FoamFile  section  are  the  matehToTolerance  and  pointSync  commands  which 
are  not  important  right  now  and  should  be  left  as  is.  Next,  you  will  see  a  patches  section,  which 
is  where  we  will  do  our  editing  to  join  surfaces  under  one  boundary  name  and  type. 

The  user  should  edit  their  createPatchDict  file  patches 0  section  to  look  like  the  screen 
capture  on  the  next  page.  The  only  patch  will  group  the  hull,  midbody ,  and  stern  surfaces  into 
one  boundary'  named  hull,  and  this  new  boundary  will  remain  a  type  wall. 

Now  we  are  ready  to  combine  the  surfaces,  but  first  it  is  generally  a  good  idea  to  copy 
our  constant  directory  before  wc  combine  our  surfaces,  so  there  is  always  a  reference.  So  for 
Linux  users,  simply: 

>>  cp  -r  constant  origeonstant 

to  copy  our  reference  constant  folder  to  orig  eonstant.  Then  in  the  case  directory  enter 
44 crcatePatch ”  in  the  command  line.  The  resulting  screen  dump  should  look  like  the  following 
screen  capture. 
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0 - 


w 

/ 

F  leld 

1 

|  OpenFOAM: 

The  Open  Source  CFD  Toolbox 

w 

/ 

0  peration 

|  Version: 

1.0 

\\  / 

\v 

A  nd 

M  anipulation 

|  Web: 

1 

http://www.openfoan.org 

FoanFile 

{ 

version 

fornat 

root 

case 

instance 

local 

c  lass 
object 


2.0: 

ascii ; 


"body  1”: 


dictionary; 
createPatcheDict ; 


II  A************************************ 

//  Tolerance  used  in  Hatching  faces.  Absolute  tolerance  is  span  of 
//  face  tines  this  factor. 
natchTolerance  IE-3; 


//  Do  a  synchronisation  of  coupled  points. 
pointSync  true; 


//  Patches  to  create. 

//  If  no  patches  does  a  coupled  point  and  face  synchronisation  anyway, 
patches 
( 

{ 

nane  hull; 
type  wall; 

construe tFron  patches; 
patches  (bow  nidbody  stern); 


): 


} 


// 


*.  +  ****  +  *****■**  +  +  **  +  **♦  +  *  +  ***♦**  +  *♦•♦+  +  ★***•*•*•*•  +  +  *♦*******♦**  +  ★  +  •*-•*•«»  +  ★  +  +  *♦**♦ 


// 
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[delaneykQamazon  bodyl_Tutorial] £ 
[delaneykQanazon  bodyl_Tutorial] £  createPatch 

/* - 


w 
\\  / 
\\  / 
w/ 


/ 


F  leld  |  OpenFOAM:  The  Open  Source  CFD  Toolbox 

0  peration  |  Version:  1.5-dev 

A  nd  |  Revision:  exported 

M  anipulation  |  Web:  http://www.OpenFOAM.org 


Exec  :  createPatch 
Date  :  Apr  08  2010 

Tine  :  15:46:37 

Host  :  anazon.dt .navy .nil 

PID  :  21526 

Case  :  /san/hone/delaneyk/NavyFOAM-1. 5-dev-rev995/delaneyk-l . 5-dev/run/bodyl/boundingI 
nProcs  :  1 


//*-***** 

Create  tine 


*  *  // 


— >  FOAM  Warning  : 

Fron  function  dlLibraryTable: :open(const  fileNaneA  funct ionLibName) 
in  file  db/dlLibraryTable/dlLibraryTable . C  at  line  86 

could  not  load  /san/hone/delaneyk/NavyFOAM-1 . 5-dev-rev995/NavyF0AM/lib/ lmux64CccDP( 
e8RASModelllprintCoef fsEv 
Reading  createPatchDict . 


Using  relative  tolerance  0.001  to  natch  up  faces  and  points 


Create  polyMesh  for  tine  =  0 

Adding  new  patch  hull  of  type  wall  as  patch  7 


Moving  faces  fron  patch  bow  to  patch  7 
Moving  faces  fron  patch  nidbody  to  patch  7 
Moving  faces  fron  patch  stern  to  patch  7 


Doing  topology  nodification  to  order  faces. 

Synchronising  points. 

Points  changed  by  average :0  nax:0 

Renoving  patches  with  no  faces  in  then. 

Renoving  enpty  patch  bow  at  position  0 
Renoving  enpty  patch  nidbody  at  position  1 
Renoving  enpty  patch  stern  at  position  2 
Renoving  patches. 

— >  FOAM  Warning  : 

Fron  function  forces :: forces(const  ob jectRegistryA  obr,  const  dictionary^  diet) 
in  file  forces/forces .C  at  line  78 
No  fvMesh  available,  deactivating. 

Writing  repatched  nesh  to  1 

End 
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You  will  notice  that  the  second  to  last  line  states  “Writing  repatched  mesh  to  I'\  The 
createPatch  command  will  write  the  new  patched  surface  information  into  a  directory  whose 
name  is  the  first  time  step  output  of  your  future  run,  which  in  this  ease  is  //. 

So  now  you  need  to  get  rid  of  the  old  constant  directory'  and  move  the  new  7/ directory  to 
constant  In  Linux  this  would  be  accomplished  by: 

»  mv  constant  orig  eonstant 

»  mv  l  constant 

Now  the  transportProperties  and  RASProperties  files  need  to  be  placed  in  the  constant 
directory.  The  transportProperties  and  RASProperties  files  must  ALWAYS  be  present  in  the 
constant  directory  when  using  simpleFoam .  In  Linux  this  would  be  accomplished  by: 

»  mv  transportProperties  constant/ 

»  mv  RASProperties  constant/ 

Open  the  constant/polyMesh/boundary  file  and  you  will  see  that  the  three  patches  are 
now  grouped  together  in  the  hull  patch. 

Now  we  have  the  geometry  imported  and  named  as  we  want  for  the  run.  A  good  next  step 
is  to  export  the  geometry'  into  a  visual  package  (EnSight,  ParaView,  etc.)  and  make  sure  that  all 
surfaces  are  grouped  and  labeled  correctly.  To  export  the  geometry,  use  foamToEnsight  for 
EnSight,  foamToVTK  for  Para  View,  and  no  additional  command  is  needed  for  ParaFoam.  So 
now  take  a  minute  or  two  and  inspect  your  geometry  in  your  package  of  choice.  Your  geometry 
should  look  like  the  Figure  below,  with  the  appropriate  surface  labels  in  the  post-processor. 


Material  Properties 

The  next  step  is  to  set  up  the  material  properties  for  the  fluid.  For  the  simpleFoam  solver, 
only  the  kinematic  viscosity  is  required  in  the  eonstant/transportProperties  dictionary  file.  Open 
the  transportProperties  file,  it  should  look  like  the  screen  capture  below.  No  editing  is  necessary, 
just  note  that  kinematic  viscosity  is  always  set  in  transportProperties. 


0* - - - *\ 


w 

/ 

F  ield 

1 

1 

OpenFOAM: 

The  Open  Source  CFD  Toolbox 

w 

/ 

0  peration 

1 

Version: 

1.3 

\\  / 
w/ 

A  nd 

M  anlpulation 

1 

1 

Web: 

http : //www . openfoar.org 

FoarFile 
{ 

version 
forrat 

root 
case 
instance 
local 

c  lass 
object 

} 

//******.*♦*..*******************♦*♦***// 

transport Model 
nu 


//  . . * . . . * . . . * . .  // 


Newtonian ; 

nu  [0  2  -1  0  0  0  0]  1.5l5le-7; 


2.0; 
ascii ; 


"constant” ; 


dictionary : 
transportPropert ies; 


The  kinematic  viscosity,  nu,  is  entered  in  as: 
nu  nu  f 0  2-1  000 J  1.5151e-7 ; 

where  [0  2  - 1  0  0  0  0]  sets  the  units  based  on  [Mass  Length  Time  Temperature  Quantity  Current 
Luminous  intensity].  The  kinematic  viscosity  dimensions  are  Length  /Time  or  in  SI  units*  m7s. 
The  value  of  the  kinematic  viscosity  is  set  to  1/Re  where  the  Reynolds  number  is  6.6M  (L  =  U  - 
1 .0).  Remember  that  this  value  must  always  be  consistent  with  the  Reynolds  number. 

Re  =  UL/nu 

0/  directory  (Initial  and  Boundary  Conditions) 

Now  we  turn  our  attention  to  the  initial  and  boundary  conditions,  which  are  stored  in  the 
0/  directory. 

For  the  icoFoam  solver  only  U  and  p  files  are  needed  in  the  0/ directory. 

Open  the  0/U  File;  it  should  look  like  the  screen  capture  on  the  follow  ing  page. 

In  the  0/U  file  notice: 

U  is  a  vector  field  quantity.  AH  U  values  must  be  set  in  vector  format,  (X  X  X). 
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The  U dimensions  must  match  the  variable  by  M,L,T,...  so  velocity  is  [0  1-1  0  0  0  OJ 

The  internalField  sets  the  initial  flow  field  condition  for  U.  For  this  case  the  fluid  is 
initially  at  free  stream  everywhere  (Ux=l  and  Uy,Uz  =  0) 

The  boundaryField  sets  velocity  boundary  conditions  for  ALL  surfaces.  All  surfaces 
must  be  included  with  proper  EC's.  All  surface  names  must  match  the 
constant/polyMesh/boundary  surface  name  exactly.  The  order  of  surfaces  is  not  important,  but 
the  names  must  match  identically. 

The  velocity  boundary  conditions  for  the  five  surfaces  arc  as  follows: 

The  hull  is  set  with 

type  fixed  Value; 
value  uniform  (0  0  0); 

to  apply  the  no-slip  boundary  condition  to  the  walls. 

The  farfield  is  set  with 

type  zeroGradient ; 

to  apply  a  zero  velocity  gradient  at  the  farfield  boundaries. 

The  inlet  is  set  with 

type  fixed  Value ; 
value  uniform  (10  0); 

for  inflow  velocity  of  Ux=\. 

The  outlet  is  set  with 

type  zero  Gradi  ent; 

to  apply  a  zero  velocity  gradient  at  the  outlet  boundary. 

The  symmetry j  is  set  with 

type  symmetry  Plane; 

All  symmetry  plane  boundary  conditions  need  to  have  type  symmetry  Plane. 

Again,  the  order  of  the  surfaces  in  the  0 /...  files  doesn’t  matter,  but  the  names  and  types 
MUST  be  consistent  with  those  listed  in  the  constant/polyMesh/boundary  file. 
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I  \\  / 

I  \\  / 

I  \\  / 

I  w/ 

V - 

FoamFile 

{ 

version 

format 

class 

location 

object 

} 

//***♦* 

dimensions 

internalField 

boundaryField 

{ 

hull 

{ 

type 

value 

} 

farf leld 

{ 

type 

} 

inlet 

{ 

type 

value 

} 

outlet 

{ 

type 

} 

symmetry 

{ 

type 

] 

} 


F  ield 
0  peration 
A  nd 

M  anipulation 


|  OpenFOAM:  The  Open  Source  CFD  Toolbox 
|  Version:  1.5-dev 
|  Web:  http ://www. OpenFOAM.org 

I 


2.0; 

ascii ; 

volVectorField ; 
"0M ; 

U; 


******************************** 

[01-10000]; 
uniform  (1  0  0) ; 


f ixedValue ; 
uniform  (0  0  0); 


zeroGradient ; 


f ixedValue ; 
uniform  (1  0  0); 


zeroGradient ; 


symmetryPlane ; 


I 

I 

I 

I 

V 


// 


// 


n 


************************************************************************* 


// 
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Now  open  the  0/p  file,  it  should  look  like  the  screen  capture  below. 

In  the  0/p  file  notice: 

p  is  a  scalar  field  quantity.  All  p  values  must  be  set  as  a  scalar,  X \  value. 

The p  dimensions  must  match  the  variable  by  M,L,T,...  so  pressure  is  [0  2  -2  0  0  0  0], 
because  in  icoFoam  p  is  actually  the  pressure  divided  by  the  density,  thus  the  SI  units  would  be 
m2/s2. 

The  internal  Field  sets  the  initial  condition  for  p.  For  this  case  (and  most  others)  we  do 
not  care  about  the  absolute  value  of  the  pressure,  p,  so  we  just  set  it  to  0  for  ease. 

The  boundaryField  sets  pressure  boundary  conditions  for  ALL  surfaces.  All  surfaces 
must  be  included  with  proper  BC’s  that  are  consistent  with  the  constant/polyMcsh/boundary 
surfaces. 

The  pressure  boundary  conditions  for  the  five  surfaces  arc  as  follows: 

The  hull  is  set  with 

type  zeroGradient; 

for  the  no  slip  wall. 

The  farfield  is  set  with 

type  zeroGradient; 

to  apply  a  zero  pressure  gradient  at  the  farfield  boundaries. 

The  inlet  is  set  with 

type  zeroGradient; 

to  apply  a  zero  pressure  gradient  at  the  inlet  boundary. 

The  outlet  is  set  with 

type  fixed  Value; 
value  uniform  0.0; 
to  set  a  reference  pressure  boundary  at  the  outlet. 

Note:  At  least  one  boundary  in  all  domains  must  have  a  set  pressure. 

The  symmetry  is  set  with 

type  symmetry  Plane; 

All  symmetry  plane  boundary  conditions  need  to  have  type  symmetry  Plane. 
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3* - *-  C++  -* - *\ 


w 

/ 

F  ield 

1 

|  OpenFOAM: 

The  Open  Source  CFD  Toolbox 

w 

/ 

0  peration 

|  Version: 

1 . 5-dev 

\\  / 

A  nd 

|  Web: 

http : //www . OpenFOAM . org 

w/ 

M  anipulation 

1 

FoanFile 

{ 


version  2.0; 

fornat  ascii; 

class  volScalarField; 

location  "0" ; 

object  p; 

) 

//*************+******+***+***********+// 


[02-20000]; 
uniform  0; 


dimensions 

internalField 

boundaryField 

{ 

hull 

{ 

type 

} 

farf leld 

{ 

type 

} 

inlet 

{ 

type 

} 

outlet 

{ 

type 

value 

} 

symmetry 

{ 

type 

} 

} 


zeroGradient ; 


zeroGradient ; 


zeroGradient ; 


f ixed Value ; 
uniform  0.0; 


symmetryPlane ; 
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Now  open  the  0/nuTilda  file,  it  should  look  like  the  screen  capture  below. 

In  the  0/nuTilda  file  notice: 

nuTilda  is  a  scalar  field  quantity.  All  nuTilda  values  must  be  set  as  a  scalar,  X,  value. 

The  nuTilda  dimensions  must  match  the  variable  by  so  viscosity  is  [0  2  -1  0  0 

0  0],  thus  the  SI  units  would  be  m2/s. 

The  interna! Field  sets  the  initial  condition  for  nuTilda .  For  this  case  it  is  set  to  le-8, 
which  is  ^10%  of  the  kinematic  viscosity. 

The  boundaryField  sets  nuTilda  conditions  for  ALL  surfaces.  All  surfaces  must  be 
included  with  proper  BC’s  that  arc  consistent  w  ith  the  constant/polyMesh/boundary  surfaces. 

The  nuTilda  boundary  conditions  for  the  five  surfaces  are  as  follows: 

The  hull  is  set  with 

type  fixed  Value; 
value  uniform  0; 

for  the  no  slip  wall. 

The  farfteld  is  set  with 

type  zeroGradient ; 

to  apply  a  zero  nuTilda  gradient  at  the  farfield  boundaries. 

The  inlet  is  set  with 

type  fixed  Value; 
value  uniform  lc-8; 

to  set  10%  of  the  kinematic  viscosity  at  the  inlet  boundary. 

The  outlet  is  set  with 
type  zero  Gradi  cn  t; 

to  apply  a  zero  nuTilda  gradient  at  the  outlet  boundary'. 

The  symmetry  is  set  with 

type  sym  m  ctryPla  ne; 

All  symmetry  plane  boundary  conditions  need  to  have  type  symmetry  Platte. 
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FoamFile 

{ 


version 

format 

class 

location 

object 


2.0; 

ascii ; 

volScalarField; 
"0" ; 

nuTilda ; 

******** 


} 

//*****♦* 

dimensions  [02-10000]; 

internalField  uniform  le-08; 


************* 


***** 


*  // 


DoundaryField 

{ 

hull 

{ 

type 

value 

} 

farf ield 

{ 

type 

} 

inlet 

( 

type 

value 

} 

outlet 

{ 

type 

} 

synnetry 

{ 

type 

) 


fixedValue; 
uniforr  0; 


zeroGradient : 


fixedValue; 
uniforn  le-08; 


zeroGradient ; 


sywietryPlane ; 


// 


************************************************************ 


************ 


// 
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Now  open  the  0/nut  file,  it  should  look  like  the  screen  capture  on  the  next  page. 

This  file  is  required  by  the  simpIeFoam  solver  when  the  Spalart-Allmaras  turbulence 
model  is  used,  but  its  turbulent  viscosity  input  values  are  not  used  in  the  calculations.  It  is  most 
likely  a  bug  in  the  code.  Nonetheless,  valid  types  and  path  names  are  required  in  the  nut  file. 

In  the  0/ nut  file  notice: 

nut  is  a  scalar  field  quantity.  All  nut  values  must  be  set  as  a  scalar,  X ,  value. 

The  nut  dimensions  must  match  the  variable  by  M,L,T,...  so  viscosity  is  [0  2  -1  0  0  0  0], 
thus  the  SI  units  would  be  m2/s. 

The  internal  Field  sets  the  initial  condition  for  nut.  For  this  case  it  is  set  to  le-6,  this 
value  is  not  important  to  the  calculation,  so  this  is  simply  a  general  ballpark  value. 

The  boundaryField  sets  nut  conditions  for  ALL  surfaces.  All  surfaces  must  be  included 
with  proper  BC’s  that  are  consistent  with  the  constant/polyMesh/boundary  surfaces. 

The  nut  boundary  conditions  for  the  five  surfaces  are  as  follows: 

The  hull  is  set  with 

type  zero  Gradient ; 

to  apply  a  zero  nut  gradient  at  the  no  slip  boundary. 

The  farfield  is  set  with 

type  zero  Gradient; 

to  apply  a  zero  nut  gradient  at  the  farfield  boundaries. 

The  inlet  is  set  with 

type  fixed  Value ; 
value  uniform  le-6; 

to  set  turbulent  viscosity'  at  the  inlet  boundary. 

The  outlet  is  set  with 
type  zeroGradient; 

to  apply  a  zero  nut  gradient  at  the  outlet  boundary. 

The  symmetry  is  set  with 

type  symm  etry Plane; 

All  symmetry  plane  boundary  conditions  need  to  have  type  symmetry  Plane. 
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3 


- *.  c++  -* - *\ 

=========  i  l 

\\  /  F  ield  |  OpenFOAM:  The  Open  Source  CFD  Toolbox  | 

\\  /  0  peratlon  |  Version:  1.5-dev  | 

\\  /  A  nd  |  Web:  http://www.OpenFOAM.org  | 

\\/  M  anipulation  |  I 


\* - 

FoamFile 

( 

version 

format 

class 

location 

object 

} 

//****♦* 

- */ 

2.0; 
ascii ; 

volScalarField; 

"O'*; 
nut ; 

*******************************  ii 

dimensions 

[02-10000]; 

internalField 

uniform  le-06; 

ooundary Field 
( 

hull 

{ 

type 

} 

farfield 

{ 

type 

} 

inlet 

{ 

type 

value 

} 

outlet 

{ 

type 

} 

symmetry 

{ 

type 

} 

} 

zeroGradient ; 

zeroGradient ; 

f ixedValue ; 
uniform  le-06; 

zeroGradient ; 

symmetry Plane; 

system/  directory  (Solver  Settings) 

Now  we  will  look  at  some  of  the  solver  settings  and  controls  that  are  located  in  the 
system/  directory.  We  will  focus  on  the  controlDict ,  fvSolution ,  fvSchemes ,  and 
decomposeParDict  files.  We  already  used  the  createPatchDict  to  merge  multiple  surfaces. 

Open  the  system/control Diet  dictionary  file.  It  should  look  like  the  screen  capture  below. 
The  controlDict  file  sets  all  of  the  run-time  parameters  and  output  information.  This  is  also 
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where  run-time  libraries  and  functions,  such  as  force  outputs  over  a  patch  and  dynamic  mesh 


libraries  are  specified. 

EF - 

- 

|  \\  /  F  ield  |  OpenFOAM:  The  Open  Source  CFD  Toolbox  | 

|  \\  /  0  peratlon  |  Version:  1.3  I 

I  \\  /  A  nd  |  Web:  http://www.openfoan.org  | 

|  \\ /  M  anipulation  |  I 

A* - V 

//  FoanX  Case  Dictionary. 

FoanFile 


{ 

version 

fornat 

2.0; 
ascii ; 

root 

case 

instance 

local 

"tutorial"; 

"bodyl" ; 

"system” ; 

c  lass 
object 

} 

die  t ionary; 
controlDict ; 

//*'**** 

libs  ( " libnavyFini teVolune . so”  " libnavylnconpressibleRASModels . so” ) 


application 

sinpleFoan ; 

startFron 

startTine; 

startTine 

0; 

stopAt 

endTine; 

endTine 

1500; 

deltaT 

1; 

writeControl 

t ineStep; 

writelnterva 1 

500; 

purgeWrite 

0; 

writeFornat 

ascii ; 

writePrecision 

6; 

writeConpression  conpressed; 


tineFornat 

general ; 

t inePrecision 

6; 

graphFornat 

raw; 

runTineModif lab le  yes: 

(Screen  capture  continued  on  next  page...  ) 
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functions 

( 

forces_Hull 

{ 

type  forces; 

//Library  to  load 

functionObjectLibs  ( “ libforces. so“ ) ; 

//Name  of  patch  to  integrate  forces  over 
patches  (  hull  ); 

//Reference  density  for  fluid  -  can  be  changed  later  ... 
rholnf  1.0; 

//Origin  for  moment  calculations 
CofR  (000); 


^  /  ★  ************************‘********ir-**********************  +  *********'*****-**  /  / 


At  the  top,  finite  volume  and  turbulenee  model  libraries  are  dynanueally  loaded  by 
libs(“..  ”); . 

The  solver  specified  in  application  input  does  not  matter.  The  solver  is  specified  on  the 
command  line  or  in  a  script  file.  Thus,  this  is  an  insignificant  line  for  our  purposes. 

The  solver  settings  are  fairly  obvious,  and  more  detail  is  provided  on  page  U-108  of  the 
User’s  Guide  (http://foam.sourceforge.net/doc/6uides-a4/UserGuide  pdf).  For  now  we  will  only 
cover  a  broad  view  of  the  file. 

We  know  that  simpleFoam  is  a  steady  solver.  Thus  the  solver  will  artificially  iterate  in 
“time”,  where  1  second  is  an  iteration.  Here  we  see  that  the  solver  start  from  startTimc  =  0,  and 
will  iterate  in  steps  of  deltaT  =  1  until  cndTimc  =  1500.  The  data  will  be  written  in  ASCII  format 
in  directories  according  to  writelntervaL  Notice  that  runTimcModifiable  is  chosen  to  yes,  this 
means  that  we  can  make  changes  to  the  control  Diet  in  the  middle  of  a  run,  and  they  will  be 
adjusted  on  the  fly,  as  opposed  to  having  the  settings  set  in  stone  for  the  whole  calculation 

One  important  note  is  that  to  start  a  calculation  from  a  previous  solution  the  startFrom 
entry  must  be  switched  to  latestTime,  and  desired  start  time  information  (directory  and  BC’s) 
must  be  present  in  the  ease  directory.  We  will  delve  into  this  further  later  on. 

Now  open  the  system/fvSolution  dictionary  file.  It  should  look  like  the  screen  capture  on 
the  next  page. 
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w 
\\  / 
\\  / 
w 


/ 


F  ield 
0  peration 
A  nd 

M  anipulation 


OpenFOAM:  The  Open  Source  CFD  Toolbox 
Version:  1.4 

Web:  http://wvfw.openfoar.org 


\* - v 


FoanFile 

{ 

version 

format 


2.0; 

ascii ; 


root 

case 

instance 

local 


class 

object 


dictionary: 

fvSolution; 


// 


// 


solvers 

{ 


p  PCG 

{ 


preconditioner 

DIC; 

tolerance 

le-7 

relTol 

0.01 

ninl ter 

1; 

maxlter 

200; 

}: 

U  PBiCC 

{ 

preconditioner  DILU; 

tolerance  le-07 ; 

relTol  0.0; 

ii  ini  ter  1; 

}: 

nuTilda  PBiCC 

{ 

preconditioner  DILU; 

tolerance  le-08; 

relTol  0.01; 

minlter  1; 

}: 


SIMPLE 

( 

nNonOrthogonalCorrectors  0; 
pRefCell  0; 
pRefValue  0; 


} 

(Screen  capture  continues  on  the  next  page...  ) 
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relaxationFactors 


{ 

p  0.3 

U  0.4 

nuTilda  0.4 

k  0.4 

omega  0.4 

} 


W  *************************************************************************  // 

The  fvSolution  file  contains  linear  solver  information  as  well  as  solver  algorithm  settings 

The  solvers  section  contains  linear  solver  settings  for  pressure,  velocity,  and  turbulent 
viscosity.  Note  that  for  this  case  we  are  using  preconditioned  conjugate  gradient  solvers  ( PCG 
for  symmetric  matrices  and  PBiCG  for  asymmetric  matrices),  but  we  also  commonly  use  multi¬ 
grid  solvers  (GAMG,  AAMG,  etc.).  The  solver  tolerance  and  relative  tolerance  settings  are  not 
important  right  now.  The  minltcr  command  sets  a  minimum  number  of  times  the  linear  solver 
will  iterate  on  a  variable.  It  is  usually  recommended  that  the  user  always  set  a  minimum  number 
of  iterations  >  0  to  prevent  the  solver  from  prematurely  not  solving  for  a  variable  (wc 
recommend  minlter  =  I ). 

Below  the  solvers  section  are  SIMPLE  algorithm  control  settings.  These  SIMPLE 
settings  are  not  particularly  useful  to  the  user  at  this  time,  so  only  a  broad  view  of  what  each 
setting  means  is  given.  Also,  note  that  the  PISO  algorithm  must  be  used  for  all  transient  solvers 
and  the  SIMPLE  algorithm  must  be  used  for  all  steady-state  solvers.  For  this  case  we  have 
nNonOrthogonalCorrcctors  set  to  0,  which  means  that  we  will  not  solve  the  pressure  equation 
more  than  once  per  iteration.  Note  for  future  runs,  if  the  pressure  residuals  arc  increasing  and  the 
solution  is  diverging/blowing  up,  nNonOrthogonalCorrectors  can  be  increased  to  iterate  the 
pressure  equation  more  and  may  lead  to  successful  solution  convergence.  Notice  that  wc  have  set 
cell  number  0  as  our  reference  cell,  where  the  reference  value  is  0.  This  is  the  reference  pressure 
for  the  incompressible  solver. 

Finally,  the  relaxationFactors  section  is  where  under-relaxation  factors  for  each  variable 
are  specified.  Typical  pressure  values  are  0.1 -0.4  and  typical  velocity  and  turbulence  quantity 
values  arc  0. 4-1.0.  Higher  values  correspond  to  quicker  solution  advancement,  but  will  be  more 
unstable  (greater  chance  of  solution  divergence). 

Now  open  the  system/fvSchemes  dictionary  file.  It  should  look  like  the  screen  capture 

below. 
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OpenFOAM: 

The  Open  Source  CFD  Toolbox 
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/ 

0 

peration 

1 

Version: 

1.3 

\\  / 

A 

nd 

1 

Web: 

http : //www. openfoan.org 

w/ 

M 

anipulation 

1 

-*/ 


FoanFile 

{ 


version 

2.0; 

fornat 

ascii ; 

root 

case 

11  11  . 

instance 

"systen** ; 

local 

class 

dictionary; 

object 

f vSchenes ; 

II  *****  * 


ddtSchenes 


***** 


************** 


// 


( 


} 


default 


gradSchenes 

{ 

default 


steadyState; 


Gauss  linear; 


} 

divSchenes 

( 

default  none; 


div(phi,U)  Gauss  linearUpwind  cellLinited  Gauss  linear  1.0; 

div(phi ,nuTilda)  Gauss  upwind; 

div((nuEff*dev(grad(U).T())))  Gauss  linear; 


laplacianSchenes 

( 

default  Gauss  linear  corrected; 

} 

interpolationSchenes 

{ 

default  linear: 

} 

snCradSchemes 

{ 

default  corrected; 

} 

f luxRequired 

( 

default  no; 

p: 

} 
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Many  of  the  fvSchemes  settings  are  not  particularly  useful  to  the  user  at  this  time,  so  only 
a  broad  view  of  the  settings  is  given  here.  For  more  detail  on  these  settings  consult  page  U- 110 
of  the  User’s  Guide. 

The  fvSchemes  file  sections  declares  the  following  settings: 
ddt  ->  time  discretization 
gradSchemes  ->  gradient  term  discretizations 
divSchemes  -)  divergence  terms  discretization 
laplacianSchemes  -)  Laplacian  terms  discretization 

interpolationSchemes  ■)  interpolation  of  values  from  cell  centers  to  cell  face  centers 
sttGradSchemes  ■)  surface  normal  gradient  evaluation  at  cell  faces 
fluxRequired  ■)  lists  fields  where  flux  is  generated  in  the  application 
Some  fvSchemes  notes: 

(1)  Because  simple  Foam  is  a  steady  solver  ddtSchemes  default  is  set  to  steadyState . 

(2)  The  div(phi,U)  term  is  the  convective  velocity  term,  and  44 Gauss  tinearUpwind 
cellLimited  Gauss  linear  7.0”  corresponds  to  2nd  order  upwind. 

(3)  The  div(phi ,  nuTilda)  term  is  the  convective  turbulent  viscosity  term,  and  “Gauss 
upwind"  corresponds  to  1st  order  upwind. 

(4)  The  div((nuEff*dev(grad(U).TO)))  term  requires  a  gradSchemes  input,  but  is  placed 
in  divSchemes.  This  is  probably  a  bug. 

For  all  of  the  fvSchemes  fields  a  default  value  can  be  specified  and  only  exceptions  to  the 
default  setting  would  need  to  be  specified,  or  default  can  be  set  to  none  which  means  that  the 
user  must  enter  all  values  for  the  appropriate  variables  themselves. 

Now  open  the  system/decomposeParDict  dictionary  file.  It  should  look  like  the  screen 
capture  on  the  next  page. 

There  is  a  lot  in  the  decomposer ar Diet  that  is  beyond  the  scope  of  this  tutorial,  but  the 
important  thing  to  notice  is  that  the  mesh  will  be  split  into  12  partitions  ( numberOJSubdomains 
12;)  using  the  metis  method. 

The  number  of  l’s  in  the  metisCoeffs  processorWeights  section  must  match  the 
number  in  numberOJSubdomains . 


109 


3 


\\  / 

F 

ield 

1 

|  OpenFOAM:  The  Open  Source  CFD  Toolbox 

\\  / 

0 

perat ion 

|  Version:  1.3 

\\  / 

A 

nd 

|  Web:  http://www.openfoan.org 

w/ 

M 

anipulation 

1 

FoamFile 

{ 

version 

fornat 

root 

case 

instance 

local 

class 

object 


//****♦ 


2.0; 
ascii ; 


dictionary ; 
deconposeParDict ; 


********** 


*  * 


// 


number Of Subdona ins  12: 


nethod  netis; 


sinpleCoef f s 

( 

n  (22  1); 

delta  0.001; 

) 


nierarchicalCoef f s 

( 

n  (11  1); 

delta  0.001; 

order  xyz ; 

) 


net isCoef fs 

{ 

processorWeights 

( 

111111111111 

); 

} 


nanualCoef f s 

{ 

dataFile 


distributed  no; 

//  *************************************************************************  if 
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Next  exeeute  the  settings  from  decomposeParDict  by  entering  4k dccomposePar M  on  the 
eommand  line. 

Upon  completion  of  the  domain  decomposition,  your  directory  will  have  twelve  new  files 
( processorO  processor  1 /),  whieh  all  correspond  to  the  decomposed  domain.  Your  ease 

directory  should  look  like  the  sereen  eapture  below. 


[delaneykQamazon  bodyl_Tutorial]$  1 
total  23 7 M 


-rw-r — r — 

1 

delaneyk 

users 

237M 

Apr 

8 

15:31  body l_Box-ASCII . fluent .cas 

drwxr-xr-x 

3 

delaneyk 

users 

21 

Apr 

8 

15:49  orig_constant 

drwxr-xr — 

2 

delaneyk 

users 

46 

Apr 

8 

16:41  0 

drwxr-xr-x 

3 

delaneyk 

users 

14 

Apr 

8 

16:45  forces_Hull 

drwxr-xr-x 

3 

delaneyk 

users 

67 

Apr 

8 

16:47  constant 

drwxr-xr-x 

2 

delaneyk 

users 

102 

Apr 

8 

16:50  system 

-rw-r— r — 

1 

delaneyk 

users 

650 

Apr 

8 

16:51  oFOAM.scp 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:52  processorO 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:52  processorl 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:53  processor2 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:53  processor3 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:53  processor4 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:53  processor5 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:53  processor6 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:53  processor7 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:53  processors 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:53  processor9 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:53  processorlO 

drwxr-xr-x 

4 

delaneyk 

users 

29 

Apr 

8 

16:53  processorll 

[delaneyk@amazon  bodyl_Tutor ial ] $ 


Running  the  Case 

We  are  now  ready  to  run  our  ease.  To  exeeute  this  ease  on  a  cluster  a  senpt  file  is  needed. 
For  example  purposes  the  script  file  oFOAM.scp  is  shown  below. 

Notiee  that  the  12  partition  mesh  will  be  run  on  3  nodes  with  4  processors  per  node.  The 
application  simpleFoant  is  also  specified  in  this  file. 

-It  is  now  time  to  run  the  job,  so  in  this  ease  we  type: 

>  >qsub  oFOA  M.  scp 

into  the  eommand  line.  A  file  named  log  will  eontain  all  of  the  run  information  that  would 
normally  be  output  in  a  sereen  dump. 

Remember  that  at  the  bottom  of  our  controlDict  file,  we  specified  a  function  named 
forcesHull  of  type  forces.  This  file  calculates  forees  over  the  pateh  specified  by  patches  ( hull 
in  our  ease),  and  plaees  them  in  a  directory  named  forces_Hull  under  a  time  File  name  that 
corresponds  to  startTime. 

Now  let  the  file  run  out  until  its  endTimc  of  1500  iterations. 


Ill 


1 


^PBS  -j  oe 

#PBS  -o  ./amazon. out 
#PBS  -e  ./amazon. err 
#PBS  -S  /bin/csh 
#PBS  -N  SA_bl 
#PBS  -1  nodes=3:ppn=4 
#PBS  -1  wallt ime=42:00:00 
#PBS  -V 

echo  "cd  to  the  directory’1 
cd  £PBS_0_W0RKDIR 

setenv  OPENFOAM.NP  12 

echo  "define  parameters  in  exec  statement" 

set  APPL I CATI0N=" simple Foam” 
set  ROOT=" . " 

set  CASE=MBody_l_Tutorial" 

echo  "The  current  shell  is  SSHELL" 

echo  "Number  of  processors:  SOPENFOANLNP" 

echo  "Executing  :  ^APPLICATION  SROOT  *CASE" 

echo  "Working  directory  :  £PBS_0_W0RKDIR" 

echo  "The  shell  limits  are:" 

limit 

echo  "Starting  executable...." 


npirun  -machinefile  i'PBS_NODEFILE  -np  $OPENFOAM_NP  ^APPLICATION  -parallel  >  ./log 

After  the  case  has  completed,  by  running  1500  iterations  open  up  the  forces  Hull/l  file. 
Let’s  just  say  for  tutorial  purposes  that  the  forces  have  not  converged  to  our  satisfaction,  and  we 
want  to  run  the  case  out  further  for  an  additional  2500  iterations. 

To  restart  the  case  make  the  following  changes  in  the  system/control  Diet  file: 

(1)  Change:  startFrom  startTime;  start  Front  latestTime; 

(2)  Change:  endTime  1500;  endTime  4000; 

-Now  restart  the  calculation  with 

»qsub  oFOAM.scp 

Notice  that  the  log  file  will  be  written  over  (so  make  a  copy  in  the  future  if  you  wish  to 
keep  the  original  log  file).  Also  notice  that  forces  are  now  being  output  under  forces/1501  file, 
and  the  original  forces  are  still  kept  under  forces/1 . 

Let  the  case  run  out  to  completion  after  4000  iterations. 

Now  open  the  log  file.  Some  observations: 

You  can  see  that  the  solver  started  from  time  equal  to  1500  and  iterated  until  4000. 

For  each  iteration  the  momentum  equation  (to,  Uy ,  and  Uz )  is  solved  first,  then  the 
continuity  equation  (p),  and  finally  the  turbulent  quantity  ( nuTilda ). 
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For  each  variable  linear  solver  we  ean  see  the  initial  residual,  final  residual,  and  the 
number  of  iterations  it  took  to  drop  from  the  initial  to  the  final  residual.  We  set  all  of  these 
tolerances  and  iteration  criteria  in  the  systcm/fvSolution  dictionary  file. 

There  are  also  continuity  error  reports. 

The  best  way  to  typically  monitor  the  solution  is  to  make  sure  that  the  velocity  magnitude 
stays  at  a  reasonable  number,  and  make  sure  that  initial  pressure  residuals  are  decreasing  or  are 
holding  steady  at  an  acceptable  value. 

The  last  line  of  the  time  iteration  produces  execution  and  clock  time  information.  This  is 
useful  in  gauging  the  efficiency  of  your  solution. 

Post-Processing 

Notice  that  there  are  many  time  directories  in  your  processor  directories.  Each  of  these 
directories  contains  output  information  for  their  respective  time  step. 

To  reconstruct  the  data  from  the  decomposed  processors  use  the  command 

>>  reconstruct  Par  -latestTime 

The  latestTime  means  only  reconstruct  the  last  time  in  the  processor*  files.  The 
command  -time  time #  will  reconstruct  for  a  specific  time  (time#)  only.  If  only  recons  tract  Par  is 
specified,  then  all  time  directories  in  the  processor *  files  will  be  reconstructed. 

To  look  at  the  post-processed  results  simply  type  the  following  commands,  depending  on 
the  post-processing  tool  of  choice: 

» foamToEnSight  -latestTime  to  look  at  the  results  in  EnSight 

>>  foamToVTK  - latestTime  to  look  at  the  results  in  ParaView 

where  the  command  -latestTime  is  used  to  only  look  at  the  results  from  the  last  output  time  step. 
To  look  at  the  results  for  all  time  steps  simply  leave  off  the  -latestTime  command,  and  to  look  at 
the  results  for  a  specific  time  (i.e.  0.005)  use  -time  0.005. 

To  look  at  the  results  in  ParaFoam,  no  additional  commands  are  needed,  simply  open 
ParaFoam  in  the  ease  directory. 

Your  results  should  look  like  the  axial  velocity  (Ux)  and  pressure  (Press)  contours  below. 
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Appendix  E:  ransFSNavyFoam  Wigley  Hull  Tutorial 

This  tutorial  involves  using  the  turbulent,  transient,  incompressible,  multi-phase  solver 
for  the  Wigley  hull.  Although  this  is  a  transient  solver,  this  case  will  NOT  be  run  time  accurate. 
Only  half  the  body  is  solved,  as  symmetry  is  assumed.  First,  we  will  go  over  pre-processing  and 
case  setup,  then  we  will  run  the  test  case,  and  finally  we  will  look  at  some  post-processed  results. 

For  more  detailed  information  on  the  OpenFOAM  code  and  settings  consult  the  User's 
Guide:  http://foam.sourceforee.net/doc/Guides-a4/UserGuide.pdf. 

Pre-Processing  and  Case  Setup 

Your  initial  directory  should  look  like  the  screen  capture  below, 
total  4. OK 

drwxr-xr-x  3  delaneyk  users  128  Jul  13  14:50  constant 

drwxr-xr-x  2  delaneyk  users  100  Jul  13  14:50  system 

drwxr-xr-x  2  delaneyk  users  72  Jul  13  14:50  0 

-rw-r — r —  1  delaneyk  users  646  Jul  13  14:50  oFOAM.scp 

constant/  directory 

In  the  previous  tutorials  the  mesh  needed  to  be  imported  into  OpenFOAM  from  a  3rd 
party  mesh  generator.  However,  for  this  case  the  mesh  has  already  been  imported,  so  you  will 
notice  the  polymesh/  folder  is  already  present  in  the  constant /  directory.  Open  up  your 
constant/polyMcsh/boundary  file,  it  should  look  like  the  following  screen  capture.  Notice  that 
the  hull  surface  is  of  type  wall  (viscous  surfaces  must  always  be  of  type  wall),  the  centcrplanc  is 
of  type  symmetryPlane  (symmetry  planes  must  always  be  of  type  symmetry  Plane),  and  the  rest 
of  the  surfaces  are  of  type  patch. 
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C  +  4 


w 
\\  / 
\\  / 
w/ 


/ 


I 

F  ield  |  OpenFOAM:  The  Open  Source  CFD  Toolbox 

0  peration  |  Version:  1.5.x 

A  nd  |  Web:  http://www.OpenFOAM.org 

M  anipulation  | 


\* - 

FoamFile 

{ 


-V 


version  2.0; 

format  ascii: 

class  polyBoundaryMesh; 

location  ’’const  ant /polyMesh” ; 

object  boundary; 


// 


*  *  // 


hull 

{ 


type 

nFaces 

startFace 

} 

centerplane 

{ 

type 

nFaces 

startFace 

} 

bottom 

{ 

type 

nFaces 

startFace 

} 

farfield 

( 

type 

nFaces 

startFace 


wall ; 
3074; 
762555; 


symmetry Plane; 
5858; 

765629; 


} 

top 

( 


type 

nFaces 

startFace 


// 


inlet 

{ 

type 

nFaces 

startFace 

) 

out  let 

{ 

type 

nFaces 

startFace 

} 

*  **************** 


patch ; 

3364; 

771487; 


patch ; 

8932; 

774851; 


patch; 

3364; 

783783; 


patch; 

2233; 

787147; 


patch; 

2233; 

789380; 


******* 


// 


Now  run  the  checkMesh  command  for  two  reasons: 

1 .  to  make  sure  the  mesh  was  imported  correctly 
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2.  to  asses  the  quality  of  the  mesh  for  the  OpenFOAM  solver 
Your  checkMesh  output  should  look  like  the  screen  captures  on  the  next  pages. 


[delaneykQarazon  wigley_tutorial]£  checkMesh 


\\  / 

F  ield 

1 

1 

OpenFOAM:  Tlie  Open  Source  CFD  Toolbox 

\\  / 

0  peration 

1 

Version:  1.5 -dev 

\\  / 

A  nd 

1 

Revision:  exported 

w/ 

M  anipulation 

1 

Web:  http://www.0penF0AM.org 

Exec  :  checkMesh 

Date  :  Jun  21  2010 

Tine  :  10:56:06 

Host  :  amazon. dt .navy .nil 

PID  :  29416 

Case  :  /san/home/delaneyk/NavyFOAM-1 . 5-dev-rev995/delaneyk-l . 5-dev/run/wigley/tutorial/wigley_tutorial 

nProcs  :  1 


//  *  *  *  *  * 

Create  tire 


*  *  // 


— >  FOAM  Warning  : 

Fror  function  dlLibraryTable : :open(const  fileNareA  functionLibNare) 
in  file  db/dlLibraryTable/dlLibraryTable.C  at  line  86 

could  not  load  /san/hore/delaneyk/NavyFOAM-1 . 5-dev-rev995/NavyFOAM/lib/linux64CccDPOpt/libnavyFiniteV 
olure.so:  undefined  syrbol:  _ZN4Foam6upwindIdE8typeNareE 
— >  FOAM  Warning  : 

From  function  dlLibraryTable : :open(const  fileNareA  functionLibNare) 
in  file  db/dlLibraryTable/dlLibraryTable .C  at  line  86 

could  not  load  /san /home/delaneyk/NavyFOAM-1 . 5-dev-rev995/NavyFOAM/lib/linux64CccDPOpt/libnavYlncorpr 
essibleRASModels . so .  undefined  symbol:  _ZN4Foaml4incompressible8RASModell lprintCoef fsEv 
Create  polyMesh  for  tire  =  constant 

T ire  =  constant 

Mesh  stats 

points:  273780 

faces:  791613 

internal  faces:  762555 

cells:  259028 

boundary  patches:  7 

point  zones:  0 

face  zones:  0 

cell  zones:  0 

Number  of  cells  of  each  type: 


hexahedra : 
prisrs : 
wedges : 
pyramids ; 
tet  wedges: 
tetrahedra : 
polyhedra : 


259028 

0 

0 

0 

0 

0 

0 


Checking  topology... 

Boundary  definition  OK. 

Point  usage  OK. 

Upper  triangular  ordering  OK. 

Face  vertices  OK. 

Number  of  regions:  1  (OK). 

Checking  patch  topology  for  multiply  connected  surfaces  ... 

Patch  Faces  Points  Surface  topology- 

hull  3074  3186  ok  (non-closed  singly  connected) 
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centerplane 
botton 
farf ield 
top 


inlet 

outlet 


5858  6105  ok  (non-closed  singly  connected) 
3364  3510  ok  (non-closed  singly  connected) 
8932  9126  ok  (non-closed  singly  connected) 
3364  3510  ok  (non-closed  singly  connected) 
2233  2340  ok  (non-closed  singly  connected) 
2233  2340  ok  (non-closed  singly  connected) 


Checking  geometry. . . 

This  is  a  3-D  nesh 

Overall  domain  bounding  box  (-4  -9.40395e-38  -4)  (12  8  1.2) 

Mesh  (non-empty)  directions  (1  1  1) 

Mesh  (non-empty,  non-wedge)  dimensions  3 

Boundary  openness  (-2 . 25653e-16  -6.20107e-16  -3. 37491e-16)  Threshold  =  le-06  OK. 

Max  cell  openness  =  3.23887e-16  OK. 

Max  aspect  ratio  =  750.787  OK. 

Minimum  face  area  =  7.99769e-06.  Maximum  face  area  =  1.44004.  Face  area  magnitudes  OK. 
Min  volume  =  3.21742e-08.  Max  volume  =  1.44.  Total  volume  =  664.872.  Cell  volumes  OK. 
Mesh  non-orthogonality  Max:  78.8601  average:  12.7201  Threshold  =  70 
Number  of  severely  non-orthogonal  faces:  18. 

Non-orthogonality  check  OK. 

<<Writing  18  non-orthogonal  faces  to  set  nonOrtho laces 
Face  pyramids  OK. 

Max  skewness  =  1.46173  OK. 


Mesh  OK. 


End 


You  will  notice  that  there  are  18  “severely  non-orthogonal  faces.”  As  has  been  mentioned 
in  previous  tutorials,  OpenFOAM’s  checkMesh  is  very  harsh.  Sometimes  it  is  not  possible  to 
ereate  a  mesh  without  any  high  aspect  ratio,  non-orthogonal,  or  skewed  cells.  In  fact,  most 
meshes  created  will  contain  bad  cells,  and  run  fine.  However,  at  some  point  (whieh  is  not 
quantitatively  clear)  the  mesh  will  be  so  poor  it  either  won’t  run,  or  it  will  take  a  long  time  to 
run.  There  aren’t  exact  guidelines  on  OpenFOAM  mesh  quality;  it  simply  takes  experience 
running  various  meshes. 

Another  good  initial  step  is  to  export  the  geometry  into  a  visual  package  (EnSight, 
Para  View,  etc.)  and  make  sure  that  all  surfaces  are  grouped  and  labeled  correctly.  To  export  the 
geometry,  use  foamToEnsight  for  EnSight,  foamToVTK  for  Para  View,  and  no  additional 
command  is  needed  for  ParaFoam.  So  now  take  a  minute  or  two  and  inspect  your  geometry  in 
your  paekage  of  choice.  Your  geometry  should  look  like  the  pictures  on  the  next  page,  with  the 
appropriate  surface  labels.  This  mesh  is  meant  for  instructional  purposes  only  as  you  will  notice 
that  the  mesh  is  very  coarse. 


1 


Entire  Domain: 


Side  View  ol  Hull; 


z 

* 


The  constant/RASProperties  file  is  the  same  as  in  the  previous  tutorials  and  will  not  be 
covered  here. 

Open  the  constant/transportProperties  File,  it  should  look  like  the  screen  capture  below. 
No  editing  is  necessary.  However,  notice  that  the  multi-phase  solver  requires  density  (rho)  and 
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kinematic  viscosity  (nu)  for  both  the  water  (phase  1)  and  the  air  (phase  2).  Additionally  the 
surface  tension  (sigma)  is  input  at  the  bottom.  The  surfaee  tension  could  probably  be  negleeted 
for  this  case  of  a  surface  ship  (sigma  =  0.0),  but  it  must  always  be  included  at  the  bottom  of  the 
transportProperties  file. 


3* - 


w 

/ 

F  leld 

1 

|  OpenFOAM; 

The  Open  Source  CFD  Toolbox 

w 

/ 

0  peration 

|  Version: 

1.3 

\\  / 
w/ 

A  nd 

M  anipulation 

|  Web; 

1 

http: //www . openf oan . org 

FoanFile 

{ 


version 

2.0; 

forwat 

ascii ; 

root 

.... . 

case 

"wigley” ; 

instance 

M  H  * 

local 

M  tl  . 

class 

dictionary; 

object 

transportProperties ; 

} 

transportModel  Newtonian; 

phasel 

{ 

transportModel  Newtonian; 

rho  rho  [1  -3  0  0  0  0  01  1000; 

nu  nu  [02-1000  0]  le-06; 

} 

phase2 

{ 

transportModel  Newtonian; 


rho  rho  [1-300000]  1; 

nu  nu  [0  2  -1  0  0  0  0]  1.48e-05; 

) 

signa  signa  [10-20000]  0.07; 

//  -*********  —  ************—******•***  —  *  —  ***************** ******  // 


The  multi-phase  solver  also  requires  a  constant/environmentalProperties  file,  which  has 
not  been  required  in  the  previous  tutorials.  This  file  contains  information  on  the  gravity  veetor  as 
can  be  seen  on  the  next  screen  capture.  An  important  note  for  free  surfaee  flow  is  to  make  sure 
that  the  gravity  and  velocity  (from  the  0/U  file)  coincide  with  the  desired  Froude  number 
according  to  the  definition: 


Fr  = 


U 

JgL 


For  this  tutorial  the  Wigley  Hull  will  be  run  at  a  Froude  number  of  0.289,  corresponding 
to  a  Reynolds  number  of  905,000. 
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\\  / 

F 

ield 

1 

|  OpenFOAM:  The  Open  Source  CFD  Toolbox 

\\  / 

0 

peration 

I  Version:  1.3 

\\  / 

A 

nd 

|  Web:  http://www.openfoaii.org 

w/ 

M 

amputation 

1 

*\ 


V 


FoanFile 

{ 


version 

2.0; 

fornat 

ascii ; 

root 

case 

“wig  ley** : 

instance 

J 

local 

M  , 

c  lass 

dictionary; 

object 

environnent a 1 Proper t ies 

} 

//**** 

6 


g  [0  1  -2  0  0  0  0]  (0  0  -9.81); 


// 


// 


// 


0/  directory  (Initial  and  Boundary  Conditions) 

Now  we  turn  our  attention  to  the  initial  and  boundary  conditions,  which  are  stored  in  the 
0/  directory.  Again,  many  of  the  basic  eoneepts  stored  in  the  0/  directory  have  been  covered  in 
previous  tutorials,  thus  only  new  concepts  will  be  covered  here.  However,  it  is  worth  repeating 
that  ALL  surface  names  in  the  0 /...  files  must  match  the  names  from  the 
constant/polyMesh/boundary  file. 

For  the  ransInterNavyFoam  solver  with  the  SST  k-omega  turbulenee  model  only  Vf  A, 
omega ,  pd,  and  gamma  files  are  needed  in  the  0/directory. 

Open  the  0/U  file;  it  should  look  like  the  sereen  eapture  on  the  next  page.  The  hull  is  set 
to  a  no-slip  boundary  condition,  the  inlet,  farfield ,  and  bottom  boundaries  are  set  to  slip 
boundary  conditions  to  simulate  the  coordinate  system  fixed  to  the  hull  as  would  be  the  case  in 
tow  tank  tests. 
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3* - *-  C++  -* . *\ 


w 

/ 

F 

ield 

|  OpenFOAM: 

The  Open  Source  CFD  Toolbox 

w 

/ 

0 

peration 

|  Version: 

1. 5-dev 

\\  / 

A 

nd 

|  Web: 

http://www.OpenFOAM.org 

w/ 

M 

anipulation 

FoamFile 


version 
format 
class 
locat ion 
object 

} 

//***** 


2.0; 

ascii ; 

volVectorField; 
"0" ; 

U; 


************ 


// 


dimensions  [01-10000]; 


internalField  uniform  (0.905  0  0); 


boundaryField 

{ 

hull 

{ 

type 

value 

} 

centerplane 

{ 

ty  v* 

) 

bottom 

( 

type 

value 

} 

farfield 

{ 

type 

value 

} 

top 

{ 

type 

) 

inlet 

{ 

type 

value 

} 

out  let 

{ 

type 

} 

} 


f ixedValue ; 
uniform  (0  0  0) ; 


symmetry Plane ; 


f ixedValue; 
uniform  (0.905  0  0); 


f ixedValue; 
uniform  (0.905  0  0); 


zeroCradient ; 


f ixedValue; 
uniform  (0.905  0  0); 


zeroCradient ; 


//  *************************************************************************  ii 


Both  the  0/k  and  0/omega  files  are  set  up  similar  to  the  previous  tutorials,  and  should 
look  like  the  following  screen  captures. 
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V 

FoamFile 

{ 

version 

format 

class 

location 

object 


dimensions  [00-10000]; 

internalField  uniform  400; 


- -  - - 

\\  / 

F  ield 

l  OpenFOAM:  The  Open  Source  CFD  Toolbox 

\\  / 

0  peration 

|  Version:  1.5- dev 

\\  / 

A  nd 

|  Web:  http://wwiv.0penF0AM.org 

w/ 

M  anipulation 

V 


2.0; 

ascii ; 

volScalarField; 
M0" ; 
omega ; 


*  *  *  * 


// 


boundaryField 

{ 

hull 

{ 

type 

} 

centerplane 

{ 

type 

} 

bottom 

{ 

type 

value 

} 

farf ield 

{ 

type 

value 

} 

top 

{ 

type 

value 

} 

inlet 

{ 

type 

value 

} 

outlet 

{ 

type 

} 

} 


zeroGradient ; 


symmetryPlane ; 


f ixedValue ; 
uniform  400; 


f ixedValue; 
uniform  400; 


f ixedValue ; 
uniform  400; 


f ixedValue ; 
uniform  400; 


zeroGradient ; 


// 


a************************************************************************ 


// 
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The  multi-phase  solver  requires  a  0/gamma  file  which  represents  the  volume  fraction 
( gamma  =  0  =  air  and  gamma  =  1  =  water).  The  solver  is  using  the  Volume  of  Fluid  (VOF) 
method  to  solve  for  both  the  air  and  water. 

Notice  that  all  of  the  bottom  and  outlet  boundaries  are  set  to  zero  gradient.  The  top  is  set 
to  inletOutlet,  which  switches  between  a  fixed  value  and  zero  gradient  condition  depending  on 
the  direction  of  flow  across  the  boundary.  The  inlet  and  farfield  are  set  to  the  calmWater 
boundary  condition,  which  keeps  the  air-water  interface  at  a  constant  height  along  the  boundary. 
The  calmWater  condition  is  especially  important  in  avoiding  artificial  waves  at  the  inlet  and  side 
of  the  domain  during  mesh  motion  calculations.  The  centerplane  is  set  to  symmetryPlane. 
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3- - ♦-  c++  -* - *\ 


I  \\  /  F  ield  |  OpenFOAM:  The  Open  Source  CFD  Toolbox  | 

|  \\  /  0  peration  |  Version:  1.5-dev  | 

I  \\  /  A  nd  |  Web:  http://www.OpenFOAM.org  | 

|  \\/  M  anipulatlon  |  I 

v - v 

FoamFile 

{ 

version  2.0; 

format  ascii; 

class  volScalarField; 

location  "0"; 

object  gamma; 

) 

//*************+**♦*****************♦*♦  // 

dimensions  [0000000]; 

internalField  uniform  0; 


boundaryField 

{ 

hull 

{ 

type  zeroCradient ; 

} 

centerplane 

{ 

type  symmetryPlane; 

} 


bottom 

{ 

type  zeroCradient; 

} 


f arf ield 

{ 

type 

valueAbove 
valueBelow 
elevation 
axis 
va  lue 

} 


calmWater ; 
0; 

1; 

0; 

z: 

uniform  0; 


top 

{ 

type 

inletValue 

value 

} 


inletOut let ; 
uniform  0; 
uniform  0; 


(screen  output  continues  on  the  next  page...  ) 
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inlet 

{ 

type 

value Above 

valueBelow 

elevation 

axis 

value 


out  let 
{ 

type 

} 

) 


calrWater ; 
0; 

1; 

0; 


z; 

uniform  0; 


zeroCradient ; 


w 


// 


Now  open  the  0/pd  file  and  notice  that  for  the  multi-phase  solver  pressure  is  in  terms  of 
the  variable  pd ,  as  opposed  to  p  for  the  single  phase  solver.  For  the  single  phase  solver  the 
pressure  (p)  is  a  relative  pressure,  whereas  a  more  precise  (pd)  pressure  is  solved  for  in 
ransFSFoam.  The  top  and  outlet  have  the  pressure  set  to  0  and  the  rest  are  at  zero  gradient  and 
symmetry  plane. 
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-  C++  - 


\\  /  F  ield 

\\  /  0  peration 

\\  /  A  nd 

\\/  M  anipulation 


OpenFOAM;  The  Open  Source  CFD  Toolbox 

Version:  1.5-dev 

Web:  http://www.OpenFOAM.org 


FoamFile 

{ 


version 

format 

class 

location 

object 


2.0; 

ascii ; 

volScalarField: 
"0" ; 
pd: 


*  +  * 


} 

//**♦++*****♦+ 

dimensions  [1  -1-20000]: 

mternalField  uniform  0; 


+  *  * 


*  *  * 


*  *  *  * 


boundaryField 

{ 

hull 

( 

type 

} 

centerplane 

( 

type 

} 

bottom 

( 

type 

} 

farf ield 

{ 

type 

} 

top 

{ 

type 
va  lue 

} 

inlet 

{ 

type 

} 

outlet 

{ 

type 

value 

} 

} 


zeroCradient ; 


symmetry Plane 


zeroCradient ; 


zeroCradient ; 


fixedValue: 
uniform  0; 


zeroCradient ; 


fixedValue ; 
uniform  0; 


// 


•\ 


V 


// 


// 


System/  Folder  (Solver  Settings) 

Now  we  will  look  at  some  of  the  solver  settings  and  controls  that  are  located  in  the 
system /  directory,  which  contains:  controlDict ,  setFieldsDict ,  fvSolution ,  fvSchemes ,  and 
decomposeParDict  files. 
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The  systcm/dccomposcParDict  dietionary  file  was  eovered  extensively  in  the 
simple  Foam  tutorial,  thus  will  not  be  discussed  here. 

Open  the  system/controlDict  dietionary  file.  It  should  look  like  the  screen  capture  below. 

The  solver  settings  are  fairly  obvious,  and  more  detail  is  provided  on  page  U-108  of  the 
User’s  Guide  (http://foam.sourceforRe.net/doc/Guides-a4/UserGuide, pdf).  For  now  we  will  only 
eover  a  broad  view  of  the  file. 

The  solver  specified  in  application  input  does  not  matter.  The  solver  is  specified  on  the 
command  line  or  in  a  script  file.  Thus,  this  is  an  insignificant  line  for  our  purposes. 

The  ransFSFoam  solver  is  a  transient  solver,  thus  it  requires  maxCo  and  maxDeltaT 
inputs  that  specify  the  maximum  possible  Courant  (CFL)  number  and  time  step,  respectively. 
When  the  adjustableTimeStcp  is  set  to  the  time  step  specified  by  deltaT  is  ignored  and  the 
time  step  size  is  ehosen  by  the  maximum  Courant  number  set  by  maxCo. 

The  maxCo  (CFL)  command  is  very  important  to  solution  stability.  There  is  no  single 
value  that  is  used  for  all  eases,  and  in  most  eases  the  user  will  start  out  with  a  low  CFL  number 
and  then  ramp  it  up  once  the  initial  solution  transients  die  out.  This  is  a  parameter  that  the  user 
will  have  to  gain  experience  over  time  to  learn  the  best  solution  strategy.  For  now  we  will  start 
with  CFL  =  5.0  and  leave  it  as  such  throughout  the  solution,  but  it  is  not  unusual  to  start  eases 
out  with  CFL  as  low  as  1  and  ramp  it  up  into  the  100’s.  The  controlDict  is  read  continuously 
during  the  calculation,  thus  the  CFL  can  be  changed  on  the  run  without  having  to  stop  the  run. 

Near  the  bottom,  finite  volume  and  turbulenee  model  libraries  are  dynamically  loaded  by 
!ibs(“.  ”);  . 

At  the  bottom  the  hullForcc  library  is  loaded,  which  calculates  forees  and  moment  over 
the  hull  surfaee.  The  moments  are  taken  about  the  point  specified  by  COR. 
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0 . - . *\ 

I  =========  I  I 

I  \\  /  F  ield  |  OpenFOAH:  The  Open  Source  CFD  Toolbox  | 

I  \\  /  0  peratlon  |  Version:  1.4  I 

I  \\  /  A  nd  I  Web:  http://www.openfoan.org  | 

|  \\/  H  anipulation  |  I 

V - V 

FoanFi le 
{ 


version 

fornat 

2.0; 
ascii ; 

root 

case 

instance 

local 

“wig  ley" : 

c  lass 
object 

} 

dictionary ; 
controlDict : 

// . 

application 

interFoan; 

startFron 

startTine ; 

startTine 

0; 

stopAt 

endTine ; 

endTine 

50: 

del taT 

0.01; 

wri teControl 

runTine ; 

wntelnterval 

10; 

purgeWrite 

0; 

writeFornat 

ascii ; 

writePrecision 

6; 

wr  iteConpression 

uncompressed; 

t ineFornat 

general ; 

tinePrecision 

6; 

runTineHodif lable 

yes ; 

adjustTineStep 

yes; 

naxCo 

S.O; 

libs  (“libnavyFiniteV'olume  .so*'  "libnavylnconpressibleRASModels.so"  "libnyDynamcFvHesh.so”) ; 
naxDeltaT  S.e-2; 

(screen  output  continues  on  the  next  page... ) 
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functions 

( 

hullForce 

( 

type  hullForce; 

//  Where  to  load  it  fror  (if  not  already  in  solver) 
functionObjectLibs  (“libhullForce. so" ) ; 


patches  (  hull  ); 

CofR  (0.5  0  0); 

) 

): 

0/ . * . * . . . . . // 


Now  open  the  system/fvSolution  dictionary  file.  It  should  look  like  the  screen  capture  on 
the  next  pages. 

The  fvSolution  file  contains  linear  solver  information  as  well  as  solver  algorithm  settings. 

Notice  that  we  arc  using  multi-grid  (GAMG)  linear  solvers  for  all  of  the  pressure  terms 
instead  of  the  conjugate  gradient  (PCG)  solvers  from  the  previous  tutorials.  The  linear  solver 
settings  and  criteria  are  explained  in  further  detail  in  the  User's  Guide.  The  important  part  of  the 
solvers  is  noticing  the  tolerance  and  relative  tolerance  (relTof)  which  determine  when  the  solver 
will  stop  iterating. 

pCorr  is  an  initial  pressure  calculation  that  is  done  before  the  first  iteration  only  for  this 
case  (as  can  be  seen  later  on  in  the  log  file).  If  the  mesh  were  moving  there  would  be  a  pCorr 
loop  for  each  iteration. 

The  transient  solver  requires  the  PISO  algorithm  as  opposed  to  the  SIMPLE  algorithm 
that  is  required  for  steady  solvers.  Most  of  the  PISO  settings  (correctors)  specify  the  number  of 
iterations  and  subiterations  for  parameters  like  velocity,  pressure,  and  gamma.  Outer  correctors 
loop  through  all  linear  solvers  (Uy  A,  omega ,  pdy  and  gamma),  non-orthogonal  correctors  loop 
through  the  pressure  equation,  and  gamma  correctors  and  subcycles  loop  through  the  volume 
fraction. 

For  future  purposes,  the  user  is  encouraged  to  change  the  various  PISO  settings  and  look 
at  the  log  file  to  sec  how  these  settings  effect  the  solution  iterations.  Solutions  on  high  quality 
meshes  will  require  less  correctors  and  subcycles,  while  for  poor  meshes  it  may  be  necessary  to 
have  more  correctors  to  achieve  a  solution.  The  higher  the  number  of  correctors  and  cycles  the 
solution  will  be  more  stable;  however,  iteration  time  will  increase  rapidly.  There  is  no  “correct" 
answer  for  each  PISO  parameter. 

The  cGamma  parameter  specifies  the  sharpness  of  the  interface  (0  -  less  sharp  and  1  = 
most  sharp).  CoGamma  refers  to  the  gamma  solution  advancement.  For  now  the  user  should 
simply  leave  cGamma  and  CoGamma  at  0  and  0.5  for  all  cases. 

For  the  time  being,  make  sure  that  your  settings  look  like  the  screen.  More  detail  on  the 
PISO  settings  is  provided  on  page  U- 117  of  the  User's  Guide 
(http  //foam. sourceforge.net/doc/Guides-a4/UserGuide. pdf). 
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OpenFOAM: 
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w 

/ 

0  peration 

1 

Version: 

1.3 

\\  / 
w/ 

A  nd 

M  anipulation 

1 

1 

Web: 

http: //www . openf oan . org 

//  FoamX  Case  Dictionary. 


FoanFile 

{ 

version 

format 

root 

case 

instance 

local 

class 

object 

) 

//•4*** 


2.0; 

ascii ; 


"wigley" : 


dictionary: 
fvSolut ion. 


solvers 

( 


pcorr  CAMC 

{ 

tolerance 

le-4 

relTol 

0; 

minlter 

1; 

maxlter 

25; 

// 


smoother 
nPreSweeps 
nFostSweeps 
nBot  tomSweeps 


DICGaussSeidel ; 

0; 

2; 

2; 


cacheAgglomerat ion  false; 
nCellsInCoarsest Level  10; 
agglomerator  f aceAreaFair ; 

mergeLevels  1; 


pd  CAMC 

{ 

tolerance 
relTol 
minlter 
maxi  ter 

smoother 

nPreSweeps 

nPostSweeps 

nFinestSweeps 


le-7 ; 

0.01; 

i; 

25: 

DICGaussSeidel ; 

2; 

2; 

2; 


cacheAgglomerat ion  false; 
nCel IsInCoarsestLevel  10; 
agglomerator  faceAreaFair ; 


mergeLevels  1; 


(screen  output  continues  on  the  next  page... ) 
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pdFinal  CAMC 
{ 


tolerance 

le-7; 

relTol 

0.01; 

mini  ter 

1; 

maxi  ter 

100; 

nVcycles 

2; 

smoother 

DICCaussSeldel 

nPreSweeps 

2; 

nPostSweeps 

2; 

nFinestSweeps 

2; 

cacheAgglomeration  false; 
nCellsInCoarsestLevel  10; 

agglomerator 

faceAreaPair ; 

mergeLevels 

): 

1; 

V  PB1CG 

{ 

precondit loner 

DI LU ; 

tolerance 

le-09 

relTol 

0.001 

mini  ter 

1; 

): 

gamma  PBiCG 

{ 

preconditioner 

DTLU; 

tolerance 

le-08 

relTol 

0; 

minlter 

1; 

): 

k  PBiCG 

{ 

precondit ioner 

DILU : 

tolerance 

le-07 

re IToI 

0.01; 

mini  ter 

1; 

): 

omega  PBiCG 

{ 

preconditioner 

DILU; 

tolerance 

le-07 

relTol 

0.01 : 

mini  ter 

1; 

}; 

PISO 

{ 


momentumPredictor 

yes 

nOuterCorrectors 

2; 

nCorrectors 

1; 

nNonOrthogona ICorrectors 

0; 
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nCannaCorr 
nCannaSubCyc les 
cGanma 
CoCanna 


relaxationFactors 

( 


pd 

0.2 

V 

0.7 

k 

0.5 

onega 

0.5 

ganna 

0.5 

l 


2: 

1; 

0; 

0.5; 


// 


Now  open  the  system/ fvSchemes  dictionary  file.  The  file  is  the  same  as  tutorials  cases 
except  for  additional  divergence  and  flux  terms.  Under  the  divSchemes  section  gamma  (phi  and 
phirb)  divergence  terms  are  needed  for  the  VOF  solution.  Also,  pd,  gamma,  and  pcorr  flux  terms 
are  required  under  the  fluxRequired  section.  Your  file  should  look  like  the  screen  capture  on  the 
next  page. 
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1  \\  /  F  ield 

|  \\  /  0  peration 

1  \\  /  A  nd 

|  \\/  M  ampulation 

\* - 

OpenFOAM:  The  Open  Source  CFD  Toolbox  I 

Version;  1.3  | 

Web:  http://www.openfoar.org  | 

1 

- v 

FoanFile 

( 

version  2.0; 

forrat  ascii: 

root  MM; 

case  "wigley" ; 

instance  "**; 

local  MM; 

class  dictionary; 

object  fvSchenes; 

} 

//**************♦*******.  **♦.******.*.*,/ 

ddtScheres 

{ 

default  Euler; 

} 

gradScheres 

{ 


default 

} 

Gauss  linear. 

divSchenes 

{ 

div(rho*phi , U) 

Gauss  linearUpwind  cellLinted  Gauss  linear  1. 

div(phi,ganna) 
div(phirb , garra ) 

Gauss  vanLeerOi; 

Gauss  interfaceCorpression; 

div(phi.k) 
div(phi ,orega ) 

) 

Gauss  upwind; 

Gauss  upwind; 

lap lacianScheres 

i 

default 

} 

Gauss  linear  corrected; 

interpolat lonScheres 
{ 

default 

} 

1  inear ; 

snGradScheres 

{ 

defaul t 

corrected; 

} 
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f luxRequired 

{ 

default  no; 

pd: 

pcorr ; 
gamma ; 

) 


g/  ****************************  ****** *4  ******* *««*«**««*««*«•  *«««««* ********  // 

Now  open  the  system/setFieldsDict  dictionary  file.  This  dictionary  can  be  used  to  initially 
set  flow  field  parameters  over  the  entire  domain.  For  now  we  will  use  it  to  set  the  domain 
volume  fraction  up  appropriately.  The  field  is  initially  set  to  air  (defaultFieldValues  setting 
gamma  0).  The  regions  section  uses  boxToCell  to  set  every  cell  within  a  box  defined  by  the 
minimum  and  maximum  rectangular  points  (which  can  extend  outside  the  domain)  to  water 
(gamma  1). 


3* - - - *\ 


\\  / 

F 

ield 

|  OpenFOAM: 

The  Open  Source  CFD  Toolbox 

\\  / 

0 

perat ion 

I  Version: 

1.3 

\\  / 

A 

nd 

I  Web: 

http:  //www .  openf oam .  org 

w/ 

M 

anipulation 

FoamFile 

{ 


version 

2.0; 

format 

as^i i ; 

root 

J 

case 

•tt*  . 

instance 

"system” ; 

local 

class 

dictionary ; 

object 

setFieldsDi 

} 

//  *  *  *  *******  *  **..♦„*♦*  *  ***.*.***♦...*♦*  n 

def aul tFieldVa lues 

( 

volScalarFieldValue  gamma  0 

): 

regions 

( 

boxToCel 1 

{ 

box  (-10  0  -20)  (20  20  0); 

f ieldValues 
( 

volScalarFieldValue  gamma  1 

): 

} 

): 


// 


// 
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The  setFields  command  will  alter  the  0/gamma  file,  so  it  is  wise  to  make  a  copy  of  your 
original  file.  This  can  be  done  in  Linux  by: 

>>  cp  —r  Of  gamma  O/Original _gamma 

Now  we  are  ready  to  set  the  initial  flow  field  with  the  setFields  command  which  executes 
the  setFieldsDict  dictionary.  Enter  setFields  into  the  command  line.  Your  output  should  look 
like  the  screen  capture  below. 


[delaneykQanazon  wigley_tutorial ] $  setFields 


\\  / 

F  ield 

1 

1 

OpenFOAM:  The  Open  Source  CFD  Toolbox 

\\  / 

0  peration 

1 

Version:  1. S-dev 

\\  / 

A  nd 

1 

Revision:  exported 

\v 

M  anipulation 

1 

Web:  http://wvw.0penF0AM.org 

Exec  :  setFields 

Date  :  Jan  21  2010 

Tine  :  11:0S:12 

Host  :  anazon.dt .navy .nil 

PID  :  31380 

Case  :  /san/hone/ delaneyk/NavyFOAM-1 . 5-dev -rev995/delaneyk-l. S-dev/run/wigley/tutorial/wigley_tutorial 
nProcs  :  1 


//****♦***♦ 

Create  tine 


// 


: — >  FOAM  Warning  : 

Fron  function  dlLibraryTable; :open(const  fileNaneA  funct lonLibNane) 
in  file  db/dlLibraryTable/dlLibraryTable.C  at  line  86 

could  not  load  /san/hone/de laneyk/NavyFOAM-1 . S-dev-rev995/Navy FOAM/ lib/ linux64CccDPOpt / 1 ibnavylnconpr 
essibleRASModels .so:  undefined  synbol:  _ZN4Foanl4inconpressible8RASModelllprintCoef f sEv 
Create  nesh  for  tine  =  0 


Reading  setFieldsDict 

Setting  field  default  values 

Setting  volScalarField  ganna 

Setting  field  region  values 

Adding  cells  with  center  within  box  (-10  0  -20)  (20  20  0) 
Setting  volScalarField  ganna 


End 

[delaneykQanazon  wigley_tutorial ] S 

It  is  wise  to  make  sure  that  the  setFields  command  did  what  it  was  supposed  to  do  by 
viewing  the  results.  Previous  tutorials  go  over  how  to  get  OpenFOAM  results  into  EnSight 
( foamToEnSight )  and  ParaView  (foamToVTK). 

Once  the  case  is  imported  into  your  post-processor  of  choice,  your  domain  should  look 
like  the  picture  below.  You  should  have  gamma  equal  to  zero  above  z=0  and  gamma  equal  to  1 
below  z=0. 
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Running  the  Case 

Now  the  problem  is  set  up  correctly  and  ready  to  run.  The  final  step  is  to  decompose  the 
domain  by  the  decomposePar  command. 

You  should  now  have  8  processor  files  (processorO  ^  processor 7)  located  in  your  case 
directory. 

The  final  step  is  to  submit  your  script  (oFoam.scp  in  this  example),  and  then  the  job  will 

run. 


>>qsub  oFOAM.scp 

Remember  that  at  the  bottom  of  our  controIDict  file,  we  specified  a  function  named 
hullForces .  This  file  calculates  forces  over  the  patch  specified  by  patches  ( hull  in  our  case),  and 
places  them  in  a  directory  named  hullForce  under  a  time  file  name  that  corresponds  to 
startTime . 

Now  let  the  file  run  out  until  its  cndTime  of  50  seconds. 

You  can  open  or  tail  the  log  file  to  monitor  your  solution  residuals  and  look  at  your  linear 
solution  strategy  that  was  set  under  the  PISO  section  in  systcm/fvSolution .  Monitoring  the 
pressure  residuals  and  the  velocity  magnitude  value  from  iteration  to  iteration  will  give  you  a 
good  idea  of  your  solution  strategy.  Rapidly  increasing  pressure  residuals  or  velocity  magnitudes 
|  over  eonseeutive  iterations  usually  mean  the  solution  is  diverging. 

Post-Processing 

Notice  that  there  are  many  time  direetones  in  your  processor  directories.  Each  of  these 
directories  contains  output  information  for  their  respective  time  step. 

To  reconstruct  the  data  from  the  decomposed  processors  use  the  command 
> >  reconstnictPar  latestTime 
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The  - latestTime  means  only  reconstruct  the  last  time  in  the  processor *  files.  The 
command  - time  time #  will  reconstruct  for  a  specific  time  (time#)  only.  If  only  reconstruct  Par  is 
specified,  then  all  time  directories  in  the  processor  *  files  w  ill  be  reconstructed. 

To  look  at  the  post-processed  results  simply  type  the  following  commands,  depending  on 
the  post-processing  tool  of  choice: 

» foamToEnSight  -latestTime  4  to  look  at  the  results  in  EnSight 

» foamToVTK  -latestTime  ■)  to  look  at  the  results  in  Para  View 

where  the  command  -latestTime  is  used  to  only  look  at  the  results  from  the  last  output  time  step. 
To  look  at  the  results  for  all  time  steps  simply  leave  off  the  -latestTime  command,  and  to  look  at 
the  results  for  a  specific  time  (ie  0.005)  use  -time  0.005. 

To  look  at  the  results  in  ParaFoam,  no  additional  commands  are  needed,  simply  open 
ParaFoam  in  the  case  directory. 

Your  results  should  look  like  the  following  figures. 

Bow  Wave: 
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Free  Surface  Contour  Plot: 
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